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SBMI- ANALYTICAL  SEMI -NUMERICAL  METHOD  TO 
TREAT  THE  THREE-DIMENSIONAL  STOKES  FLOW 


/1 40 


Lin  Shengtian  Vu  Wangyi 


(Department  of  Mechanics,  Beijing  University) 
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SUMMARY  This  article  makes  use  of  the  axially  symmetrical 
Sampson  fldw  body  internal  distribution  method  to  numerically  solve 
three  dimensional  S.tokes  flow  problems.  As  far  as  the  complicated  . 
integrals  which  appear  in  ..this  method  are  concerned *  we  initially 
obtained  systematic  and  simple  results  of  succesive  estimation  .type 
analysis.  Because  of  this*  the  high  accuracy  of  the  method  is 
guaranteed*  and  the  method  has  the  advantage  -of  a  small  amount  of 
required  calculation.  Ve  carried  out  specific  calculations  for  the 
three  dimensional  Stokes  flow  lateral  prolate  spheroid  and  Cassini 
oval  flow  bodies  in  a  uniform  oncoming  flow.  Ve  obtained  pressure 
distributions  for  object  surfaces  and  drag  coefficient  values  with 
very  good  convergence  and  accuracy.  As  jin  actual  example  of  the 
Stokes  flow  for  a  three  dimensional  body*  this  article  alBO  calculated 
the  Stokes  flow  for  a  three  dimensional  prolate  ellipsoid  in  a  uniform 
oncoming  flow.  Moreover*  it  obtained  for  the  first  time  the 
corresponding  drag  coefficients  and  pressure  distributions. 

Creep  flow  for  viscous*  non-compressible  flow  bodies  has  always 
been  a  problem  having  basic  significance  and  engaging  people's 
interest.  Following  along  with  the  products  in  the  last  twenty  years 
of  such  fields  as  biological  engineering  and  chemical  engineering  as 
well  as  with  the  development  of  the  fields  of  science  and  technology* 
the  problem  of  the  Stokes  flows  around  the  viscous  flow  bodies 
surrounding  any  configuration  of  object  has  undergone  rapid 
development.  People  have  hoped  to  find  a  type  of  method  which  had 
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relatively  broad  applicability ,  as  veil  as  a  relatively  high  degree  of 
accuracy  and  a  reasonable  amount  of  calculation  involved.  Such  people 
as  Gluckman  [1]  selected  for  use  a  ring-form  singularity  point  surface 
distribution  method  to  handle  unbounded  Stokes  flow  movements  for  any 
convex  shaped  axially  symmetrical  body.  Youngren  and  Acrivos  1.2)  have 
brought  forward  the  method  of  continuously  distributing  the  Stokes 
flow  body  on  the  surface  of  the  object  in  order  to  handle  unbounded 
Stokes  flow  movements  for  three  dimensional  bodies.  Both  these  two 

types  of  methods  include  numerical  value  integr  ations.  Because  of 
this,  the  amount  of  calculation  is  relatively  large.  Most  recently, 
such  people  as  Wu  Wangyi  f3,4,5,6,7)  have  presented  a  new  method  which 
takes  Sampson  flow  bodies  and  distributes  them  continuously  on  a 
certain  selected  line  segment  of  the  body  internal  axis  of  symmetry 
and  on  a  certain  selected  central  surface  of  the  body  interior.  This 
successfully  solved  the  problem  of  the  Stokes  winding  flow  around  any 
body  symmetrical  around  its  long  axis  and  around  any  body  symmetrical 
around  its  side  axis.  This  method  achieved  thoroughly  ideal  results. 
Ve  calculated  the  axially  symmetrical  flow  movements  for  such  bodies 
as  spheroids  ,  flattened  spheres,  elongated  Cassini  oval  bodies, 
flattened  Cassini  oval  bodies,  and  long  bodies  with  pointed  tips.  It 
goes  without  saying  that  these  reflected  the  special  characteristics 
of  the  bodies  as  a  whole  in  terms  of  drag  coefficients  or  localized 
pressure  distributions  which  appear.  All  of  these  quantities 
possessed  relatively  good  convergence  results  and  relatively  high 


calculation  accuracy.  Therefore,  this  explains  the  fact  that  Sampson 
flow  bodies  present  a  type  of  fast,  highly  accurate  numerical  value 
method  appropriate  to  axial  points  of  singularity  in  order  to  solve 
for  the  Stokes  flow  movements  for  any  configuration  of  axially 
symmetrical  body 

The  test  diagrams  in  this  article  generalize  a  method  suggested 
in  reference  (3)  for  handling  problems  of  three  dimensional  Stokes 
winding  flows.  We  found  the  algebraic  expressions  for  several 
collections  of  functions  and  several  successive  estimation  equations 
in  equations  expressing  the  Sampson  progression,  and  completed  several 
totally  different  integrals  in  the  coordinate  directions.  They  were 
expressed  in  the  form  of  analytic  successive  estimation  formulae. 
Therefore,  the  success  of  the  method  was  guaranteed.  For  the  sake  of 
the  effectiveness  of  inspection  methods,  we  did  research  into  the 
three  dimensional  Stokes  problem  of  a  spheroid  in  a  uniform  oncoming 
flow  from  the  side  direction.  This  used  very  little  in  the  way  of 
position  points  and  a  small  amount  of  calculation.  We  obtained  quite 
good  drag  coefficients  and  pressure  distribution  values  in  agreement 
with  precise  solutions.  This,  therefore,  confirmed  the  convergence 
characteristics  and  precision  of  this  method  in  solving  three 
dimensional  problems.  After  this,  we  selected  for  use  the  same  type 
of  method  and  calculated  the  drag  coefficient  value  for  the  lateral 
direction  three  dimensional  winding  flow  for  the  Cassini  oval  body  of 
rotation  to  act  as  a  practical  example  of  the  treatment  of  the  concave 
cavity  portion  of  three  dimensional  flow  movements.  Finally,  we  did 
research  on  the  Stokes  flow  movements  of  three  dimensional  ellipsoids, 
obtaining,  for  the  first  time,  drag  coefficient  values  and  pressure 
distribution  values  for  flows  around  three  dimensional  ellipsoid 
bodies  in  a  uniform  oncoming  flow  from  any  given  direction. 
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I.  Lateral  Direction  Stokes  Winding  Flows  Around  Bodies  Symmetrical 
About  the  Long  Axis 
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Consider  the  Stokes  flow  (See  Fig.  1)  around  the  -v  axis  of  a 
body'of  rotation  in  a  uniform  oncoming  flow  along  the  direction  of  the  ’ 
axis  at  a  velocity  U  .  Select  L%U*  \MlLt  to  respectively  be  the 
characteristic  parameters  for  length,  velocity  and  pressure.  In  this,  L 
is  a  certain  characteristic  length  of  the  body,  n-  is  cue  dynamic 
viscosity  coefficient  of  the  flow  body.  Copying,  the  method  used  in 
reference  (3),  on  the  line  connecting  the. rate  of  curvature  centers  A 
and  B,  we,  separately  or  in  a  continuous  manner,  position  the  Sampson 
flows  along  the  direction  of  the  *  axis. 

1  .  The  Separated  Line  Distribution  Method 


S 


On  the  basis  of  certain  considerations,  on  the  line  segment  AB, 
we  select  ^  different  points.  Their  coordinates  are  (0,  £i>,  /« i»  •••» 

.  At  these  points,  we  position  Sampson  flows.  After 

the  disturbance  fields  produced  by  them  are  added  to  the  uniform 

oncoming  flow,  it  gives  us 
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U(r,  cc.,$.(r,  £,)+D.,T.Cr,  I,)]  _ 


(1  .1) 


In  this 


U=(v.-ltv„v„  p- pj  =  (ritz) ,  F  l'\rlt  z)*Rk,Fl"  -J-  ,o) 


3  p  c  i  > 


O',.*)) 


In  this,.  r  is  a. vector  radius.  ri  “V  (x~i, )•'+!/• +z‘  ,  /?,=>/(*-£,)* +y*  »• 
v„v„v,  are  respectively  the 

projections  of  the  velocity  vectors  on  the  x,y»*  axes.  P~  is 
the  pressure.  P  is  the  pressure  at  an  infinitely  distant  point.  C.„  D„ 
are  parameters  awaiting  determination.  The  equations 
expressing  F.“*(,=- i,2,s,6)  arej 


F.,1,<r,?)=r-‘«-‘,Cp.C£)+2 /.(C)] 
/;,.‘0<rf<f)«(ii+i)r-»  /.M(C) 


O',*)  **  (»  +  l)r- «•-»  , (C) 


(1  .2) 


-2*  r- 


/.(C) 


At  these  positions  /?-✓  *'+(/'  ,  r  =>/  ,  C  =  */r,  .  /, 

are  respectively  "  order  Legendre  multi-term  equations  and  n| 
order-1 /2  degree  Gegenbauer  multi-term  equations. 
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Take  a  cross  section  at  the  quantity  jV-m  in  the  infinite 
progression  in  equation  (1.1).  After  that,  select  points  AfAr  on  uhe 
surface  of  the  object.  Make  the  viscous  adhesion  conditions  at  the 
points  A/ A'  i  satisfactory.  Obtain  the  equation  for  ?.  MV  -  and 
determine  the  2  unknown  parameters  and  and  C,K  D„ 

After  this,  substitute  into  equation  (1.1)  and  obtain  an  equation 
approximately  expressing  velocity  and  pressure.  The  drag  coefficient 
is 


S  Dtl 

*  "  i'.s LU  0*3) 


2.  Continuous  Line  Distribution  Method 

On  AH  ,  we  continuously  distribute  the  Sampson  flow  in  the 
direction  of  tha  s  axis.  We  obtain  the”  equation  expressing  the  flow 
field,  which  is: 


v  •  ^ 

■  z1*?  J,/[C.(!)$.<r,*)  +£>„(!)T,<r,!} 


]</£ 


(1.4) 


In  this,  CM)*  DM)  are  strength  distribution  fractions  awaiting 
determinate  ion.  £  at  poinis  A ,  B  'equals  d  ancl  -d  .  ■•...■  .r.u 

divide  it  into  A/  sections.  In  each  section,  obc  o t.  z..0^n 
distribution  function  can  be  approximated  by  the  use  ox  an  n  order 
multi-term  equation  <»“0,  l,  •••)  .  On  the  basis  of  the  experiences 

in  reference  151*  the  accuracy  of  approximations  by  second  order, 
multi-term  equations  is  high,  and  the  amount  of  calculation  involved 
is  reasonable.  Because  of  this,  this  article  mainly  makes  use  of 
parabolic  approximation.  Due  to  this  fact,  in  each  segment,  we  select 
the  two  end  points  and  the  midpoint  to  be  the  interpolation  points. 

The  strength  distribution  in  each  section  is  approximately  replaced  by 
the  use  of  the  second  degree  functions  at  each  extrapolation  point. 
Moreover,  we  take  the  cross  section  of  the  infinite  series  at  the 
quantity  iV  +  i  .  Due  to  equation  (1.4),  it  is  possible  to  write: 
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(1  .5) 


CC.,d)  S,  (r,  $>  +  D.(  (S)T.(r,4)]</£ 


In  this 


(1  .6) 


Moreover, 

AT,,-  -  j—a-dli)U-d,i) 

>  ^<~d,l~dii,  HI~diJ  —  </|(,  dit 


Here,  </,,,  c/j.  are  respectively  the  coordinates  for  the  two  end 
points  and  the  midpoint  of  each  small  section. 

We  take  equation  (1.6)  and  substitute  it  into  equation  (1.5). 
Moreover,  we  make  the  integral  transformation  .  After  going 

through  the  manipulations,  we  obtain 


.V  ♦  t  If 

2  S  2  )!).„ ,k p,., * + 

•  - 1  ( •  t  » - 1 


(1.7) 


In  this 


c\vu,  6\'; o> 
9.,. -(o.-i  ..  esr,.., 
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/ 143 


/ ?l»>  _ 

“•l  l.»“ 

2  9i.*C /.*.';**  ix-d,,y,x)^/il 

t  •  1 

1 

”<i  (.*■ 

2  9i.tU U\+,(x-d,  ,y,z) - 

/:r 

<•1 

i  % ; 

G « « ».  — 
“ • .  i.j" 

v  2  9‘.hi/i!'il,(.x-d,  ,y,z) 
!•« 

_  _  <*-<*, .X^t-d^Xx-d,,)  , 

*“ - - “A» 

9t.t^i(*-^u)-~C(*-dtl)  +  (x~dtt)  +  (x—dll)])  •  A* 
«».*■  A* 


A»« 


UThTHT 


Moreover , 


/•<*> 


</x  i-1,2  A-0,1.2," 

dx  i“5,  (  n«l,2,3,  ••• 


(1  .8) 


The  success  or  failure  of  the  current  method  is  determined  by  the 
ability  or  inability  to  complete  the  integration  in  equation  (1.8). 

It  is  particularly  fortunate  that  we  successfully  went  through  a 
series  of  iterative  approximation  formulas  and  solved  for  an  equation 
expressing  the  closed  form  of  all  integrals  in  (1.8). 

First,  we  make  use  of  the  related  equations 
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Sc 


P  <l>  _  *  P  ID  _  H~  1  1  P  (  I  > 

TV  7T  *•"* 


xt**1**0^""  TrF"i]*'  ^  -£-</?' +*')F.‘'i 


Fi*1-  — Fill 

II 


F;»«  -A. 


F**'  =  </?'  +  *l>FJ”- 

VI 


(1.9) 


We  obtain 
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<ft>  _ 
*•  t 


4ft) 

•  -  l.ft 


*-0,1.2.-  (1.10) 
n«i,2,S*.. 


-  -"-ZA  /-.' ‘ r.1  (t/,  +  r*)/'!‘. 


The  problem  reduces  to  one  of  solving  for  an  equation  to  express 
It  is  not  difficult  to  derive  the  successive  approximation  /; J\ 
formula  shown  below: 


/ 


U> 

•t  l - 


>  > 


.< k-JHn-2 ) 

n(n-l)  ""  1  5 


**2,3, 

#=2, 3, ... 


(1  .11  ) 


If  one  is  able  to  solve  for 
then,  from  equation  (1.11), 
easy  to  demonstrate  that 


f 

f 


«.  ») 

<  i  v  , 

•  *  * 
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Wv  > 

•  t  J  » 


V'-t  i  ■  ‘  /  -V)  • 

/  I  ,  5  »  /  0,  5  » 


is  entirely  determined. 


It  is 


10 


•)  h— 1,2,  • 

VI 


(1.12) 


b) 

f  «•>  „ 

Hziy'+’z’) 

C<«- 

r'lMy’+SsM/i!*,.,- 

3(»— 

2)* 

+  («-3)/:i 

1,$^ 

xRFil\l 

fir,- 

xz 

/<•»«  1  .  JL 

/  «•> 

_  '  XZ 

y  6 

yr  * 

/ 1.  > 

riy'+s1) 

O 

/.i. » 

*»-» 

A-.,  :  . 

*  • 

-i- /.r,~  J-,„(£±i.) 


0.13) 


(1  .14) 


0  .15) 


In  c)  and  d) 


I?-* 


Prom  the  successive  approximation  formula  given  below  and  the 
equations  expressing  B„  B ,  ,  one  determines  that 


£«  =  In(rH-x),  F,  =  r 


0  .16) 


Arriving  at  this,  and  with  the  help  of  equations  (1 . 1 0) — ( 1 .1 5) ,  it  is 
very  easy  to  derive  all  the  integrals  in  equation  (1.8). 


In  the  equation  expressing  the  physical  quantities  in  formula 
(1.7),  there  are'  2lV(2A/  +  l>  •  unknown  functions.  In  orcter  to  solve 

for  them,  let  the  velocity  field  at  these  points  satisfy  the 
non-slippage  movement  condition.  Then,  obtain  the  linear  algebraic 
equation  set  of  2N12M  + 1)..  order,  which  gives  precise  values  for  C.i  »  D»> 
We  make  use  of  the  method  of  solving  for  the  inverse  of  Cm„  D„ 
any  type  of  matrix  and  solve  for  the  numerical  values  in  this  set  of 
equations.  We,  then,,  obtain  C.„  D,j0  .  After  substituting  into 
equation  (1.7),  we  then  obtain  equations  which  approximately  express 
velocity  and  pressure. 

The  formula  for  the  drag  coefficient  is  /145 


2  :  "  .  •  »  ■. 

n’f  ff  S  +4  .,<4-1  di ) 

• i - 1  4 


(1  .17) 


II.  Lateral  StokeB  Plows  Around  Spheroids 

In  order  to  check  the  convergence  characteristics  and  accuracy  of 
the  method  discussed  above,  we  first  consider  the  unbounded  Stokes 
flow  movements  in  a  lateral  direction  around  the  .*  e^ia  of  spheroid 
bodies  of  rotation  and  carry  out  a  comparison  of  precise  solutions. 

In  a  rectangular  coordinate  system,  the  equation  for  spheroidal 
object  surfaces  is 

x^ttcosAt  y  mb  sin  A  tin  0,  z**bsiuAcosO 


a^b,  (See  Pig.  l) 

In  these  problems,  the  precise  solutions  for  object  surface 
pressure  distributions  and  drag  coefficients  are  respectively  [81 


A  - 
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c'  2e+  (3e’  —  1)  In 
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i-#J 
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p-p,.. 

fiUla 
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In  this,  emy/a -~b‘  is  focal  distance,  and  #»c/a  is 

the  eccentricity. 

Vhen  the  separated  line  distribution  method  is  chosen  for  points 
of  singularity,  we  select  the  long  axis  .*  as  the  characteristic 
length.  (_c»  ®»  0)  are  the  points  B  .  Due  to  the 

fact  that  flows  around  the  body  are  completely  symmetrical  about  the 
plane  :*-o  ,  the  singularity  point  strength  relating  to  the  plane  *= 

is  aldo  symmetrical.  On  the  x  axis,  we  select  for  use  the 
bisection  principle  to  position  M  points  of  singularity.  Vhen  .tf»i 
,  a  point  of  singularity  is  positioned  at  the  origin  point,  '/hen 
M>2  »  the  position  of  the  point  of  singularity  on  the  x  axis  is 


\fC-i  O-D  <1-1.  M), 


Vhen  the  continuous  line  distribution  method  is  selected  for  use 
with  the  points  of  singularity,  we  take  M  and  divide  it  into  c 
equal  sections.  Moreover,  we  choose  the  interpolation  value  points  to 
be  the  center  points  of  the  various  sections.  Because  of  this,  we 
have 


V  i  O'-ih  V  d, 


"  2 


</'-!.  "’iV), 


For  the  placement  of  points  on  the  surface  of  spheroids,  we  make 
use  of  the  polar  angle,  equal  division  principle.  On  the  boundary 
line,  A-O-xlz,  .  This  possesses  special 

importance.  The  reason  for  this  is  that  it  is  precisely  these  which 
determine  the  spheroid  projection  area  perpendicular  to  the  direction 
of  movement.  Vhen  we  do  a  detailed  examination  of  the  parameter 
matrix  of  the  set  of  linear  equations,  we  discover  that  the  points  of 
singularity  are  on  these  points.  In  order  to  overcome  the 
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difficulties  discussed  above,  we  use  four  nearby  points  (*/ 

(xl2)-6t  .  which  respectively  take  the  place  of  0  =  O,*/2;l=*Q,jr/2. 

Here  d„d„  slt  are  constants  awaiting 

definition,  Their  -selection  must  follow  this  type  of  principle  *  that 
is,  that  it  is  both  necessary  to  avoid  the  influence  of  singularity 
and  also  necessary,  as  much  as  is  possible,  to  control  the  windward 
surface  of  the  object.  In  the  calculations  in  this  article,  we  choose 
0.05  degrees  or  0.1  degrees. 

The  results  for  drag  coefficients  which  are  obtained  through  the 
use  of  the  calculation  methods  discussed  above  are  set  out  in  Table  1 
and  Table  2.  Moreover,  at  the  same  time,  these  tables  present  a 
comparison  of  precise  solutions  for  various  data.  It  is  possible  to 
see  that,  in  the  few  circumstances  in  which  -V/»  X  do  not  match  well 
with  the  position  points,  it  is  possible  to  obtain  drag  coefficient 
values  which  are  quite  satisfactory  and  which  have  good  convergence 
characteristics  and  accuracy.  In  the  dispersed  singularity  point  /1 46 

distribution  method,  when  the  length  to  width  ratio  6/a<0.7  ,  one 

convergence  characteristics  and  accuracy  of  drag  coefficients  are 
relatively  good.  It  Is  possible  to  achieve  four  significant  digits. 
However,  when  the  width  to  length  ratio  is  larger  than  0.7,  the 
convergence  characteristics  of  drag  coefficients  turn  bad.  After  the 
strength  distribution  function  selects  for  use  a  division  segment 
which  approaches  a  parabola,  the  difficulties  above  are  thoroughly 
resolved.  In  calculating  the  width  to  length  ratio  straight  through 
from  0.9  to  0.01 ,  the  convergence  characteristics  and  the  accuracy  of 
the  results  get  better  and  better.  Estimates  for  width  to  length 
ratios  somewhat  smaller  than  this  show  that  this  type  of  good 
convergence  result  will  continue  straight  through  as  the  process  is 
continued. 

Table  1  and  Table  3  present  the  pressure  distributions  on  the 
surface  of  spheroids  as  calculated  using  all  the  methods  discussed 
above.  Moreover,  at  the  same  time,  it  presents  a  comparison  of  the 
precise  values  calculated  for  these  materials.  These  results 
demonstrate  that  the  results  of  calculations  made  using  the  parabolic 
section  approximation  method  also  possess  relatively  good  accuracy  as 
concerns  such  localized  characteristic  quantities  as  pressure  dist¬ 
ribution.  Approximate  and  precise  calculations  for  pressure  agree 
with  each  other  to  4  or  5  significant  digits. 
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Table  1  The  Separated  Line  Distribution  Method  1 .  Spheroid  Drag 
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Table  3  Pressure  Distribution  on  Spheroid  Surfaces  (  Continuous  Line 
Distribution  Method)  1  .  Precise  Solution 

The  sample  calculations  above  explain  the  fact  that  the  parabolic  /i  47 
line  segment  approximation  method  offers  a  type  of  method  for  handling 
three  dimensional  Stokes  flow  problems  which  not  only  has  the 
advantages  of  fast  convergence,  high  accuracy,  and  savings  on  time, 
but  is  also  not  limited  by  width  to  length  ratiios.  In  the  continuous 
line  distribution  method  calculations  below,  parabolic  segment 
approximation  is  always  used. 
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Table  4*  Cassini  Oval  Body  1.  Drag  Coefficient  2.  Surface  Pressure 
Distribution 
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Table  5  Three  Dimensional  Ellipsoid  Drag  Cofficients  1  .  * MY  is  the 
point  separation  number  in  the  direction  of  separation  on  xhe  central 
surface. 


III.  The  Three  Dimensional  Stokes  Flovs  of  Rotating  Cassini  Oval 
Bodies 

In  the  section  above,  ve  vent  through  the  Stokes  flov  for  a 
spheroid  in  a  lateral  uniform  oncoming  flov  and  looked  at  the 
convergence  characteristics  and  the  accuracy  of  three  dimensional  flov 
movements  as  determined  by  the  Sampson  singularity  point  distribution 
method.  The  results  vere  entirely  satisfactory.  This  section  vill  go 
a  step  further  in  selecting  the  methods  discussed  above  in  order  to 
determine  the  Stokes  flov  in  a  uniform  oncoming  flov  from  the  z 
direction  around  any  given  form  of  rotation  body  in  the  *•  direction. 
Ve  choose  the  Cassini  oval  body  of  rotation  to  be  an  example.  When 
the  parameters  are  not  the  same,  the  Cassini  oval  body  presents  many 
types  of  forms.  There  are  completely  convex  ones.  There  are 
partially  convex  ones.  And,  there  are  partially  concave  ones.  This 
is  very  useful  for  the  increased  applicability  of  inspection  methods. 

The  equation  for  a  Cassini  body  of  rotation  in  the  *•  direction, 
in  a  three  dimensional  rectangular  coordinate  system,  is 
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*“peo*^,  v—p  tin  4  *»psin^co*0 

P—  c*  cos  2  #+ ✓o*-c*  sin1 24  («>c>0) 


In  this,  p,  4  are  respectively  vector  radii  and  polar  angles  in  a  /148  • 
polar  coordinate  system.  c  is  the  focal  length.  a  v;  some 
constant.  ,6»v'o,-c*  is  the  coordinate  tor  one  point  of 

intersection  between  the  oval  body  and  the  ft!  (ose  Pig.  2). 

Talc?  6 -to  be  the  characteristic  length.  Diagrams  of  the  meridian 
surfaces  of  the  Cassini  oval  bodies  when  e-*o.8;  2.5,  so  are  shown 
in  Pig.  3.  c»0.8  corresponds  to  a  fully  convex  form.  c~2,5 

a  slightly  concave  section.  c~so  has  a  very  deep  concave  dip 
section.  Obviously,  the  larger  is  the  more  severe  is  the  concave 
dip.  We  first  select  A«p.«/p.i.  -  %/i+2c’  to  be  the  parameter. 

Then,  we  consider  the  several  situations  where  U5,  1.3,  2.0,  7.0  e 

Hs  position  a  point  of  singularity  at  (~c,c).  f  and 
respectively  make  use  of  the  separation  line  distribution  method  and 
the  parabola  line  segment  approximation  to  carry  out  calculations. 

The  singularity  point  placement  and  the  handling  of  the  position 
points  are  the  same  as  In  the  case  of  spheroids. 


Pig.  2 
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Pig.  3 

IV.  Stokes  Plows  Around  Three  Dimensional  Ellipsoids  in  Uniform 
Oncoming  Plows 


Concerning  the  Stokes  flows  around  an  ellipsoid  in  a  uniform 
oncoming  flow  from  a  certain  direction,  it  is  possible  to  solve 
respectively  for  the  oncoming  flow  around  the  ellipsoid  in  the 
directions  of  the  three  coordinate  axes  .  Not  to  lose  universality, 
it  is  possible  only  to  consider  the  flow  around  an  ellipsoid  in  a 
uniform  oncoming  flow  in  the  r  direction.  We  position  a  Sampson  fl^w 
on  the  central  plane  in  the  body  perpendicular  to  the  oncoming  flow 
(take  the  focal  length  to  be  half  the  axis  of  the  ellipse).  Make  the 
points  of  singularity  in  the  *  direction  (or  V  direction)  distribute 
continuously.  Moreover,  we  use  separate  distribution  in  the 
direction.  Because  of  this,  it  is  possible  to  make  use  of  formulae 
(1.1)  and  (1 . 5 ) — ( 1 .6)  from  the  first  section.  The  position  points  on 
the  surface  of  objects  all  make  use  of  the  principle  of  the  equal 
division  of  polar  angles.  In  making  use  of  the  methods  described  in 
the  two  sections  above,  we  avoided  the  singularity  of  the  linear 
equation  set  parameter  matrix. 


/1 49 


The  results  obtained  from  calculations  for  drag  coefficients  are 
presented  in  Table  5«  The  calculation  results  show  that,  within  a 
certain  range  of  width  to  length  ratios,  segments  in  a  certain 
direction  approximate  parabolas.  In  another  direction,  we  selected 
for  use  the  separation  approximation  Sampson  singularity  point  surface 
distribution  method.  This  was  sufficient  to  obtain  drag  coefficient 
values  of  over  three  significant  digits.  This,  therefore,  explains 
the  applicability  of  the  Sampson  singularity  point  iterative  addition 
method  in  handling  problems  involving  th*»  Stokes  flow3  around  three 
dimensional  bodies. 
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Abstract 


•  ^  „  j  •  t  •  I  i  •  •  i*  *•*£!  *  ^ 

This  paper  deals  with  the  numerical  solutions'  of  the  three-di 

.•r M 


uuension 


Stokes  flow  by  means  of  the  axisymijiet^ic  Sampsorileis  ^listribujipn  method 
whithin  the  body.  The  systematical  and  simple  formulas  arc  obtained  for 
the  first  time  to  evaluate  the  complex  integrals  emerged  in  this  paper. 
This  leads  to  the  advantages  of  proposed  method)  high  accuracy  and 

•  •  ■  • ./  ,i  >*■  .1  ■  •  ■  * 

small  amount  of  the  computational  work  The-thrcu-dimenjion  Stokes  flow 

!  .  .  :i  . :  ■  •  ».i.  *«  *-  • .  ,<t  .;j  , 

for  the  prolate  spheroid  and  the  revolutional  Cassini  oval  passed  by  uni- 

'  .  ......  .  .  .  ..T'S  : 

form  flow«in  the  side  direction  are.  calculated.  1  he  good  convergence  prop¬ 
erties  tin:  high  accuracy  of  thu  results  for  the  drag  factor  and  the  pres¬ 
sure  distribution  are  demonstrated.  Finally  in  this  paper  the  three-dimen¬ 


sion  Stokes  flow  passed  an  cllipsoid^arc' considered  and  the  corresponding 

drag  factor  and  the  pressure  distribution;  aao-  precepted.  . 

■  *  »  ■  «.  <  -‘i..  -**  f  [-  .  f. 
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SUMMARY  This  article  develops  a  type  of  systematic,  linear 
numerical  value  conformal  method,  which  makes  it  possible  to  take  the 
exterior  area  of  a  given  curved,  multi-sided  form  containing  several 
angular  points  and  do  a  linear,  conformal  mapping  of  it  as  the 
exterior  area  of  a  unit  circle.  Due  to  the  fact  that,  in  this  method, 
use  is  made  of  the  FFT  technique,  computational  speed  is  very  fast. 
This  type  of  method  can  be  used  in  calculating  two-dimensional, 
noncoifipressible  flows,  around  a  given  body,  as  well  as  in  figuring  the 
apparent  mass  coefficient  of  two-dimensional  bodies,  producing  the 
orthogonal,  conformal,  body-attached  coordinate  grid  required  in 
limited  difference  calculations.  Results  with  calculations  from 
several  types  of  exterior  forms  show  that,  on  the  grid  surface  pro 
duced,  the  machine  time  required  is  an  order  of  magnitude  smaller  than 
that  for  the  Thompson  method,  and  the  orthoganality  of  the  grid  is 
preserved  everywhere.  Besides  this,  in  quite  a  few  problems,  the 
governing  equation  on  the  transformation- plane  is  particularly  simple. 

I.  Introduction 

In  the  last  ten  years,  the  technology  of  producing  coordinate 
grids  that  adhere  to  surfaces  has  rapidly  developed.  This  has 
overcome  the  difficulties  of  handling  curve  boundary  conditions  in 
limited  difference  methods.  The  use  of  equations  for  elliptical  forms 
to  produce  grids,  as  developed  by  Thompson  II],  is  relatively  good  and 
flexible.  However,  the  expenditures  of  time  involved  are  relatively 
great.  If  this  method  is  used  in  calculations  for  supersonic  flow 

speeds,  the  time  occupied  in  the  production  of  the  grid  and  the  time 
used  in  the  calculation  of  the  flow  field  can  be  of  the  same  order  of 
magnitude.  Object-adhering  coordinate  grids  produced  for  conformal 
mapping  have  several  advantages.  They  preserve  orthogonality.  In 
potential  flow  problems  outside  of  shock  wive  areas,  in  areas  where 
one  often  get3  abrupt  changes  in  flow  movements  (sucn  as  the 
vicinities  of  ’'he  fo- -rac’d  and  rear  edges  of  airfoils)  one  obtains 
denser  grid  patterns.  Besides  this,  the  governing  equations  for  many 
problems,  in  the  conformal  mapping  calculation  planes,  have 
particularly  simple  forms.  Because  of  this,  it  is  widely  used  for 
calculations  in  transonic  potential  flows  [2,31  and  other  flow 
movement  problems. 
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Conformal  mapping  methods  were  used  very  early  to  handle 
non-compressible  flows  around  given  airfoils  [41 .  In  the  1970's  FFT 
[13l  was  developed  as  the  basic  method  for  fast,  conformal  mapping 
f5,6j,  which  causes  an  order  of  magnitude  increase  in  the  speed  of 
calculations.  If  the  exterior  form  included  several  points  of 
angularity,  we  generally  chose  for  use  a  series  analytic 
transformation  to  eliminate  the  angular  points  [7,3 1.  After  this,  we 
carried  out  a  transformation  of  the  smooth  boundary.  Recently,  Davis 
(9),  also  taking  Schwarz-Christoff el  transformations  as  a  basis 
developed  a  type  of  capability  to  handle  the  conformal  mapping  of 
curved  edge  exterior  forms  containing  angular  points.  The  conformal 
mapping  which  is  used  by  Jameson  [21  also  takes  FFT  techniques  as  a 
basis.  In  this,  it  is  possible  to  take  any  airfoil  and  map  it  as  a 
unit  circle.  The  rear  edge  of  the  airfoil  is  equivalent  to  an  angular 
point.  This  article  will  take  the  latter  type  of  method  and 
generalize  it  to  include  any  curved,  multi-sided  form  containing  many 
points  of  angularity.  The  method  developed  is  not  only  capable  of 
being  used  to  produce  coordinate  grids  adhering  to  bodies-  It  is  also 
capable  of  being  used  to  calculate  incompressible  flows.  It  can  also 
be  used  to  calculate  apparent  mass  coefficients.  It  is  the  foundation 
for  calculating  the  aerodynamic  derivatives  for  guided  missiles  and 
aircraft  [10-12]. 


II.  Conformal  Mapping  Methods 

The  problem  is  to  take  the  exterior  of  a  certain  curved, 
multi-sided  form  (Fig.  1)  on  physical  plane  Z_  y.nd  map  ic  oonfromaily 
as  the  exterior  of  a  unit  circle  on  complex  plane  <?i.  We  assume  that 
the  transformation  derivative  can  be  expressed  as  /ipl 


Pig.  1  Physical  Plane  z  and  Transformation  Plane  «  I .  Plane 


In  the  equation,  the  large  brackets  enclose  the  Schwarz-Christoffei 
factor.  This  is  used  to  cause  the  rate  of  slope  of  the  lines  of 
intersection  to  be  discontinuous  at  the  points  of  angularity.  M  is 
the  number  of  points  of  angularity.  e,  is  the  included  angle  at  the  /. 
th  point  of  angularity.  <?'■* •  is  the  corresponding  complex 
coordinate  for  points  of  angularity  on  the  °  plane.  They  are  on  the 
unit  circle.  is  unknown  before  the  fact.  c.=»a.-i6,  is  an 

unknown  complex  parameter.  The  number  of  quantities  in  the  expansion 
formula  N  must  be  chosen  to  be  adequately  large. 

Take  equation  (l)  and  expand  it  in  an  infinitely  distance 
neighborhood.  In  order  to  make  the  transformation  have  a  single  value 
(  that  is,  the  closed  exterior  form  mapping  being  a  closed  ext  eri.o> 
form),  L >1  the  expansion  of  the  equation,  it  is  necessary  not  to 
include  the  quant ity  l/tf  .  Because  of  this,  the  value  of  parameter  C( 
must  necessarily  be 


(2) 
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Making  use  of  this  equation,  the  expansion  formula  can  be  written  a3 


.3) 


In  order  to  determine  the  complex  parameter  -c«»  ,  on  tu-=  unit 
circle  »  a~ e"'*»  •  Because  of  this,  da^-iv  •* 

Correspondingly,  dZ’*  e'*  ds  a  .  Here,  .  **;  is  the  arc  length 
calculated  in  a  counterclockwise  direction.  ft  is  one  angle  of  slope 
of  the  tangent  or  intercept  line  (Pig.  1).  Vie  taxe  these  equations 
and  substitiute  them  into  formula  (1).  We  take  the  logarithm  and 
separate  the  real  and  imaginary  sections,  obtaining  a  precise 
definition  of  tk»  formulas  for  a«  nl  b. 
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(5) 


In  this,  Hi.0,-0)  i3  the  unit  step  function.  When  t the 

value  *■  ;•  v  •  iji  O.  Whan  _  take  the  value  to  h  .*  1  . 

As  far  as  the  exterior  :o.*  i  3.1  m  n  on  the  2  pi  n  ;  is  concerned,  ft 

acts  as  a  function  of  r  .  /J(s)  is  already  known.  v3  take  the 

circumference  of  the  unit  circle  on  the  ,T  plane  (  <0=0  to  2.*)  ; 

and  divide  it  into  J  equal  sections.  The  iterative  substitution 

procedure  for  determining  the  parameters  "»  .vud  >  ■  •  qiven 

below: 

1  .  A33US9  that  the  point  divi  1  value  s  >  •;  already 

known;  that  is,  assume  the  functional  raia-ti lirtshij?  .  s(l})\  For 
example,  in  the  initial  approximation,  it  is  possible  to  assume  this 
to  be  a  linear  relationship.  Por  the  function  S(0y,  ■  we  use 

interpolation  valuer  for  values  of  0  ietermined  for  points  of 
angularity,  that  is,  'ot,  j-  l.  2 ,—,M„ 
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2.  Prom  the  already  known  values  of  (t).  and  •<<)).  ,  the  left 

side  of  equation  (5)  becomes  the  already  <mown  function  of  i  .  It  is 

possible  to  expand  it  as  a  Fourier  series,  specifying  a,,  '■■■,  o.v  >nd 
*  »  h  •  Due  to  the  fact  that  the  exterior  form  has 

points  of  angularity,  the  function  ,  therefore,  is 

discontinuous  at  the  points  of  angularity.  It  is  possible  to 

demonstrate  that  the  second  term  of  the  left  side  of  equation  (5)  will 

produce  a  corresponding  discontinuity  at  the  points  of  angulariity. 

This  causes  the  difference  in  the  second  term  of  the  left  side  to  be  2* 
the  periodic  discontinuity  function  of  1  &  .  Because  of  this,  the 

right  side  of  this  type  of  equation  is  a  Fourier  series  with  good 
convergence. 

3. We  take  the  values  determined  for  a.  and  6 ■  and  substitute  them 
into  the  right  side  of  equation  (4).  In  this,  <»«  and  ;.re 
calculated  according  to  formula  (2).  Reconciling  this  to  the  Fourier 
series,  we  get 


(6) 


In  this,  the  value  of  /(#)  at  the  equally  spaced  poinh  0,  on  the 
circumference  of  the  unit  circle  is  already  known.  We  again  make  use 
of  numerical  value  integration  methods  to  obtain  a  new  function 
relationship  where  $  functions  as  e 


$—  «•» 


(7) 


Due  to  the  fact  th-- +,  the  overall  length  of  the  circumference  of  the 
exterior  form,  s*  - ,  on  the  plane  Z.  is  already  known,  it  is, 
therefore,  necessary 


SH---- 


C*o 


'*  f{0)d0 


(3) 
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Prom  this  equation,  we  make  a  precise  solution  for  a,»  . 

It  is  necessary  to  point  out  that,  when  making;  a  numerical  value 
integration  of  formula  (7),  if  the  domain  of  infestation  ^*<3 

includes  i  •  !tself  a  point  of  angularity  9,  ,  ,;aen,  0,<0,<0.n«, 

Acc  :  •  i  '  to  equation  (.4) »  the  singularity  of  f(0)  in  the  vicinity 
Of  i.3  .  In  this  -«</*.  •••'hen  .T<e,<2» 

,  although  i'  n  infinitely  large  singularity,  it  is  still 
possible  to  integrate  .it.  Now,  let  us  assume  that,  in  the  small 
segement  f(9)  any  possibly  approximate  the 

expression 


/<0)  =  /<0;)<0,-0>,7(0^-0^>• 


Using  an  analytic  method  of  integration,  it  is  possible  to  find 

/<0)d0~/<0,>  <0,-0,  )/(i-  •£-) 

The  same  type  of  method  leads  to  the  integral  of  the  area  (0„  0j.,). 

Because  of  this,  /153 

[*'"  n&)dO~  (*’  +  (*•♦'  =  /(f)«rOt/(t,)(ll,rtl) 

*-  "  C2^-)  '>■ 

This  equation  makes  clear  the  fact  that,  when  e, -- 2.t  ,  it  cannot 

be  integrated.  Because  of  this,  this  method  is  nob  able  to  handle 
this  type  of  situation. 

If  CO.,  0,  ,1  does  not  contain  a  point  of  angularity  or 

for  a  point  of  angularity,  then,  using  a  simple,  ladder  type 
formula  we  obtain 


I*'*'  /<0>d0-  ^-CM> +/<*»,, 
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4.  Take  the  newly  determined  s(0)  and,  using  a  low  relaxation, 
orthogonal  relationship  sK$i,  assumed  in  the  first  step,  repeat  the 
first  to  the  third  steps  in  cycles  of  iterative  substitution. 

Continue  this  until,  the  difference  between  i<0>  as  determined  in  two 
neighboring  iterative  substitutions  is  within  the  nominal  range  for 
permissible  differences,  in  the  process  described  above,  it  is  nec¬ 
essary  several  times  to  take  the  already  known  function  and  expand  it 
to  a  Fourier  series  as  well  as  taking  the  already  known  :  Fourier 
series  and  reconciling  them.  FFT  techniques  are  very  effective 
techniques  when  used  in  carrying  out  these  operations  Up].  When 
using  FFT,  the  point  separation  or  scanning  number  .  /  -just  be  chosen 
as  a  whole  number  power  of  2,  that  is  24~16,  2  -32,  2*=G4,  2’*=i28  , 

and  so  on.  The  N  chosen  in  equation  (1)  is  N-JI 2.  .  The 

example  calculation  which  was  done  demonstrates  that  the  iterative 
convergence  discussed  above  is  fast.  In  the  case  of  airfoils,  6-7 
iterative  substitutions  is  sufficient.  Jn  relatively  difficult 
problems,  the  relaxation  factor  in  step  4  can  be  chosen  as  0.5  or 
smaller. 


III.  Orthogonal  Conformal  Grid  Creation  Methods 

In  the  previous  sections,  we  have  already  taken  the  mapping  of 
the  exterior  portion  of  objects  on  the  i>lane  to  be  exterior 
portions  of  the  unit  circle  on  plane  a  .  Due  to  the  fact  that  the 
exterior  domain  is  infinitely  large,  we  use  the  transformation  re¬ 

take  the  mapping  of  the  exterior  portion  of  the  unit 
circle  plane  o!  to  be  within  the  unit  circle  plane  r  .  At  this 
time,  an  infinitely  distant  point  corresponds  to  the  center  of  the 
circle.  Within  the  unit  circle  plane  r  ,  in  the  direction  r  ,  r. ; 

make  use  of  the  interval  as  Ar  (it  is  also  possible  for  the 
intervals  to  be  unequal).  Concentric  circles  at  this  interval  as  well 
as  lines  radiating  in  the  *  direction  at  interval  -\0  creates  an 
orthogonal  grid  (Fig.  2).  This  type  of  grid  takes  the  mapping  to  be 
an  orthogonal,  conformal  coordinate  grid  adhering  to  the  object  on  its 
surface  plane.  In  order  to  determine  the  grid  points  of  the 
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Pig.  2  Diagram  of  Grid  Formation  1  .  Plane  2.  th  Line 


coordinates,  we  first  determine  the  transformation  modulus  H*  on  t’->e  r 
plane, 


LI  =  !  dZ  dZ  I  ji  dej,  -'j  i  -dZ  I 

1  «fr  i  (jv.  I  [  c(r  j  -r*  UdcfcL:.. 


(9) 


In  formula  (1),  let  .  Take  the  logarithm.  Then,  if  we 

take  its  real  component,  we  obtain  /154 


la*  ^ 

,0‘  7? 


2  x(l” -J-)***  «l“reo.(«-W  +  r» 

+  >;  rVo.  cot  h fl  + 6.  sin  n fl) 


>} 


(10) 


The  triangular  series  on  the  right  side  of  the  equation  still  uses  FFT 
technology  for  reconciliation  or  summation.  From  this,  we  find  grid 
point  positions  I dZ/da  and  H  values. 
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this  is  dus  to  the  faot  that  the  transformation  modulus  H  is  the  ratio 
between  the  grid  step  length  at  a  location  on  the  2  pi*.ne  and  the 
grid  step  length  of  the  corresponding  point  on  che  r  plane.  If  we 
have  already  found,- on  the  r-’  plane,  the  *  tu  concentric  circle  r~kSr  t 
the  corresponding  curve  on  the  physical  plane  £_ks  well  as 
the  positions  of  the  grid  points  on  it,  then,  due  to  the  orthogonal 
grid,  it  is  possible,  from  the  grid  points,  to  draw  a  perpendicular 
line  toward  the  outside  (for  example,  point  P  }  (Fig.  2).  If  we  take 
PR—H&r,  then,  it  is  possible  to  obtain  the  grid  point  on  a 
line  defined  by  *  (  (r~(A+i)Ar)  )  below.  Because  of  this,  the 

orthogonal  conformal  grid  on  the  2  plane  is  able  to  be  gradually 
created  toward  the  outside  along  a  perpendicular  direction  beginning 
from  the  surface  of  the  plane.  It  is  necessary  to  point  out  that, 
during  the  creation  process  discussed  above,  the  radial  direction  step 
length  .Sr  -lust  not  be  chosen  too  large. 


IV.  Calculation  of  Apparent  Mass  Parameter 

In  order  to  cause  the  derivative  of  the  mapping  equation  to  be 
one  at  an  infinitely  distant  point,  we  introduce  C~exp(c„)cr,  .id 

substitute  it  into  equation  (2).  After  integrating,  we  obtain 


In  this 


I,  —  gL1: 


(12) 


The  other  parameters  can  also  be  calculated,  however,  i  c 

is  very  difficult  to  express  them  as  simple  analytic  equations.  Thr 
exterior  form  mapping  of  the  exterior  plane  is  a  circle  in  plane  C  . 

Its  radius  is  J?,«exp(o«>.  .  According  to  the  definition 

apparent  mass  coefficient  in  reference  [101  as  well  as  its  symbols  aid 
calculation  formulas,  the  calculation  formulas  for  the  translation  of 
apparent  mass  coefficients  in  the  direction  of  the  horizontal  axis  .-.s 
well  as  in  the  direction  of  the  vertical  axis  are 
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««,  -  1*  /»|  -R*(At)  J 

■Ki  ■  mtl  m  -  t*/>  /««(^,) 

m  2*  p  [#}  -  A-  +  R«A,)  ] 


In  these  equations,  *  I.-  flow  body  density.  S.  i*  the  area  of  the 
exterior  form  o?  the  2  pl.^ne.  R *  .-a\  lm  -.-xH  the  symbols  for  tne 
real  and  imaginary  p^rts.  Calculations  of  the  apparent  mass 
parameters  as  well  as  *■«  corresponding  to  rotations  are 

relatively  difficult.  It  is  necessary  to  make  full  use  of  the 
parameters  developed  in  formula  (11).  /155 

V.  Sample  Calculation 


Example  1 :  Semicircle  with  radius  1 .  Pig.  3  is  fi  following 
changes  in  the  arc  length  *  .  During  the  calculations,  the  scanning 

or  spacing  number  of  the  unit  circle  on  the  plane  r  is  taken  *s  64(A0 -2*/64> 
The  radial  direction  scanning  number  is  taken  as  MKArao.i) 

In  this  example,  the  iterative  substitution  relaxation 
factor  is  taken  as  0.5  (if  it  is  taken  as  0.8,  one  does  not  get 
convergence.)  After  13  rounds  of  iterative  substitution,  the  accuracy 
obtained  is  close  to  that  for  two  iterative  substitutions,  and  the 
maximum  deviation  of  the  scanning  arc  length  is  smaller  than  0.0005. 

A  comparison  of  the  calculated  results  for  apparent  mass  and  precise 
solutions  appears  in  Table  I.  Pig.  4  is  the  orthogonal,  conformal 
grid  which  is  produced  (what  is  shown  in  the  Pig.  is  based  on  only  ;2 
points  mapped  per  cycle).  Pig.  5  is  the  pressure  sUrengtn  paramt t*=r  C, 
distribution  for  noncompressible,  horizontal  winding  flows  and  agrees 
with  the  results  of  Davis  f9l. 
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Pig.  3  Changes  in  the  Semicircular  Shape,  Tangential  Slope  Angle  a 
According  to  Arc  Length 

Pig.  4  Orthogonal,  Conformal  Grid  of  Semicircular  Outer  Section  Area 


Pig.  5  Semicircular  Shape  Noncompressible  Potential  Plow  Surface 
Pressure  Strength  Distribution 


Example  2:  Equilateral  triangle  of  height  1  .  For  the 
calculation,  choose  A0»2x/64,  Ar*0.l»  •  Choose  the  iterative 

substitution  relaxation  factor  as  0.3*  15  iterative  substitutions 

reaches  the  mandated  accuracy.  The  results  of  calculations  of 
apparent  mass  are  seen  in  Table  I.  In  the  case  of  an  equilateral 
triangle,  the  apparent  mass  coefficients  for  any  given  direction  are 
all  equal.  The  minute  difference  between  and  m*»  calculated  in 
this  example  is  caused  by  calculation  error.  Fig.  6  is  the  orthogonal 
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0.1161 

0.1453 

0 

0 

Table  I  1.  Semicircle  of  Unit  Radius  2.  This  Article  3-  Precise 
Solution  (Reference  [ll])  4.  Right  Triangle  of  Height  1  5-  Precise 

Solution  (Reference  [10,11]) 


conformal  grid  produced.  Fig.  7  is  the  pressure  strength  distribution 
for  noncompressible  horizontal  winding  flows. 

Example  3:  A  wing-fuselage  composite  assembly  cross  section 
surface.  The  half  wing  deployment  equals  2.  The  round  fuselage 
radius  is  1  .  This  example  contains  6  points  of  angularity.  Two 
interior  angles  are  zero.  Four  interior  angles  are  3.r/2,  .  For 

these  calculations,  choose  y-i28i  •  Choose  the  relaxation  factor 
as  0.15*  Fig.  8  is  a  comparison  of  the  calculated  pressure  strength 
distribution  for  horizontal  noncompressible  winding  flows  and  the 
results  of  precise  solutions.  From  this  figure  one  can  see  that,  when 
equidistant  points  on  the  unit  circle  of  plane  a  are  projected  on  the 


physical  plane,  in  the  vicinity  of  points  of  angularity  at  which  the 
interior  angle  is  greater  than  n  ,  the  points  are  usually  very  thin, 
influencing  the  accuracy  of  calculations. 


Fig.  6  Orthogonal,  Conforms!  Grid  for  Bquilateral  Triangle  Fo.’m 


Around  Equilateral  Triangle  Shape 

Pig.  8  Noncompressible  Potential  Plow  Pressure  Strength  Distribution 
Around  Wing-Puselage  Composite  Body  Cross  Section  1 .  This  Article  2. 
Precise  Solution  3.  Fuselage  Surface  4.  Wing  Surface  5.  Arc  Length 
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Example  4:  GA  ( W) — 1  supercritical  airfoil.  We  have  already  done 
specialized  work  in  calculating  the  low  speed  pressure  strength 
distribution  for  a  given  airfoil  and  have  composed  a  calculation 
program.  These  airfoils  contain  only  one  point  of  angularity.  In  /1 57 

equation  (l),  if  c,  is  taken  to  be  an  appropriate  value  (see  reference 
[2]),  it  is  possible  to  take  the  rear  edge  non-closed  airfoil  mapping 
to  be  a  closed  unit  circle.  Pig.  9  is  a  comparison  between  the 
surface  pressure  strength  distribution  calculated  when  the  angle  into 
the  wind  is  a -4.17°  and  experimental  values  for  it.  When  doing 
calculations,  one  circumference  of  the  airfoil  is  taken  as  128  points, 
that  is,  Atf =2.7/128  .  6  iterative  substitutions  reach  the 

required  accuracy.  The  whole  process,  on  the  IBM  4541  machine,  takes 
only  three  seconds  of  CPU  time.  Discrepancies  between  calculations 
and  experimental  results  are  mainly  due  to  the  fact  that,  in 
calculations,  there  is  no  consideration  of  the  effects  of  viscous 
boundary  layers.  We  have  also  compiled  a  program  to  produce 
orthogonal,  conformal  coordinate  grids  which  conform  to  objects  such 
as  a  given  airfoil.  On  the  IBM  4541  machine,  producing  a  64x16  grid 
only  required  three  seconds  of  CPU  time. 


Pig.  9  GA  (W)-1  Airfoil  Low  Speed  Pressure  Strength  Distribution 
1. Upper  Surface  2.  Lower  Surface  5-  Experiment  4-  This  Article's 
Calculations 
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A  NEW  TECHNIQUE  OF  NUMERICAL  CONFORMAL 
MAPPING  AND  ITS  APPLICATIONS 

Huang  Mingke 
•T>' an  ling  Aeronautical  Institute ) 

Abstract 

A  systematic  method  of  numerical  conformal  mapping  is  developed  in 
this  paper  to  map  the  exterior  of  an  arbitrary  2-1)  body  having  several 
discontinuities  of  surface  slop  onto  the  eiterior  of  a  unit  circle.  The  use 
of  the  1  f  T  technique  makes  the  execution  of  the  computation  very  fast. 
I  he  method  developed  can  he  applied  to  the  compulation  of  incompressible 
flow  past  an  arbitrary  2-D  body,  to  that  of  the  apparent  mass  coeffi¬ 
cients  of  the  cross-sections,  and  to  the  generation  of  the  orthogonal,  con¬ 
formal  grid  used  in  finite  difference  method.  Several  examples  presented 
show  that  the  present  method  has  an  order  faster  in  computation  of  the 
grid  grr.ei  ar  it n  than  that  of  the  usual  Thompson’s  method,  and  that  the 
grid  generated  has  the  orthogonality  everywhere.  The  another  superiority 
is  that,  for  many  problems,  the  governing  equation  usually  takes  the 
especially  simple  form  on  the  mapping  plane. 
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NUMERICAL  SOLUTION  OP  TRANSONIC  SMALL 
DISTURBANCE  PRESSURE  EQUATIONS  USING 
A  MIXED  DIFFERENCE  METHOD 


Wang  Lixia  Luo  Shijun 


(Northwestern  Polytechnical  University) 


SUMMARY  This  article  presents  a  numerical  value  method  for 
directly  solving  the  transonic  steady  state  small  perturbation 
pressure  equation 

^  l-  Mi  -  Mi  «  )««+«„»- 

As  far  as  research  into  interference  problems  with  various  tunnel 
walls  are  concerned,  compared  to  traditional  velocity  potential 
equations,  the  use  of  pressure  equations  to  be  the  governing  equations 
for  solving  transonic  flow  fields,  with  the  boundary  conditions  of 
Dirichlet  form,  makes  it  easier  to  handle  these  problems.  Moreover, 
taking  pressure  as  the  variable  to  be  solved  for,  it  is  possible  to 
reduce  accumulated  error  and  increase  the  accuracy  of  the 
calculations. 

This-  article  chooses  for  use  a  mixed  difference  method  18]  to 
solve  pressure  equations  through  numerical  value  testing  methods, 
solving  the  appropriate  difference  forms,  and  iterative  substitution 
linearization  methods.  Convergence  solutions  with  these  methods,  when 
compared  to  the  solutions  obtained  by  taking  velocity  potential 
equations  as  the  governing  equations,  give  results  which  compare 
relatively  well.  This,  therefore,  demonstrates  the  feasibility  of  the 
methods  in  this  article 

Finally,  this  article  presents  typical  sample  calculations  for 
applying  the  methods  of  this  article  to  calculating  tunnel  wall 
interference  during  transonic  wind  tunnel  tests  of  airfoils  and  the 
calculating  of  pressure  distribution  for  the  exterior  forms  of  certain 
airfoils . 


In  recent  years,  the  theory  of  self-correcting  wind  tunnels  and 
various  types  of  methods  for  the  correction  of  interference,  which 
have  been  developed,  present  quite  an  optimistic  outlook  for  the 
solving  of  tunnel  wall  interference  problems.  Among  these,  there  are 
.quite  a  few  methods  which  require  the  use  of  pressure  distributions 
which  are  measured  on  the  surface  of  wings  and  in  the  vicinity  of 
tunnel  walls  in  order  to  calculate  the  boundary  conditions  for  flow 
fields  inside  wind  tunnels  (l-7J.  This  type  of  problem  usually  only 
involves  model  surfaces  and  flow  field  pressure  distributions.  The 

methods  which  have  been  used  in  the  past  (2,5-7),  always  start  by 
taking  the  empirically  measured  pressure  boundary  condtions  and  using 
them  for  potential  velocity.  After  that,  they  did  numerical  value 
solutions  for  velocity  potential  equations.  Then,  they  turned  back 
around  and  took  the  solutions  for  velocity  potential  and  turned  them 
into  pressures.  In  this  way,  on  the  one  hand,  they  handled  the 
complexity  of  the  procedures  for  boundary  conditions,  and,  on  the 
other,  in  the  two  stage  pressure-velocity  potential-pressure 
transformation  process,  they  inevitablly  added  cumulative  error.  With 
reference  to  the  particular  characteristics  of  these  types  of 
problems,  if  it  is  possible  to  directly  solve  equations  using  pressure 
as  a  variable,  changing  the  boundary  conditions  to  Dirichlet  boundary 
conditions,  then,  it  is  possible  to  overcome  the  drawbacks  discussed 
above.  Because  of  this,  this  article,  in  dealing  with  small  transonic 
perturbations,  derives  small  perturbation  pressure  equations.  These 
equations,  in  transonic  ranges,  have  the  same  forms  as  are  found  in 
small  perturbation  velocity  potential  equations.  They  can  also  have 
mixed  forms.  These  two,  in  their  mathematical  natures,  have  certain 
points  in  common.  However,  because  of  the  appearance  of  the  term  .•<;  , 
the  pressure  equations  are  slightly  more  complicated  than  the 
velocity  potential  equations.  Moreover,  the  conversion 
characteristics  and  stability  characteristics  of  numerical  value 
solutions  to  pressure  equations  are  also  more  difficult  to  rigorously 
demonstrate  mathematically.  Because  of  this,  this  article  makes  use 
of  numerical  value  tests  to  research  the  actual  calculation  methods 
used  in  mixed  difference  methods  f 8 ]  for  solving  pressure  equations, 
and  we  have  initially  obtained  satisfactory  results 
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II.  Transonic  Steady  State  Small  Perturbation  Pressure  Equations  /160 


The  binary  transonic  small  perturbation  velocity  potential 
equation  is: 


In  this,  Af„  represents  the  rnach  number  of  pertubation  gas 

flows  which  have  not  yet  passed  through.  stands  for  the  velocity 

of  the  perturbation  gas  flows  before  they  pass  through.  ^  stands  for 

the  small  perturbation  velocity  potential.  V  stands  for  the  heat 

adiabatic  index  of  the  atmosphere. 

Take  (1)  and  solve  for  *  and  let  .  Then,  the  equation 

which  satisfies  the  partial  perturbation  velocitywin  the  x  direction 

is: 


( i -mi-  mi u ) y+*.  ml u;  -0  ( 2) 

If  we  take  note  of  the  relationship  between  the  small 
perturbation  pressure  parameter  c,  and  u  : 


C,=  - 


2  u 

K' 


(3) 


then,  as  soon  as  we  know  u  ,  we  can  conversely  solve  for  CM  . 
Because  of  this,  this  article  takes  equation  (2)  and  calls  it  the 
transonic  small  perturbation  pressure  equation. 

III.  Boundary  Conditions 

1  .  Boundary  Conditions  for  the  Calculation  of  Plow  Fields  Inside 
Wind  Tunnels 
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Pig.  1  Boundary  Conditions  for  Calculations  of  Plow  Fields  Inside 
Wind  Tunnels  1.  Airfoil  2.  Side  Pressure  Control  Surface  3*  Gas 
Permeable  Wall 


As  Pig.  1  shows,  when  calculating  flow  fields  inside  wind 
tunnels,  it  is  possible  to  select  the  calculated  domain  to  be  the 
rectangular  area  shown  in  the  Pig.  The  boundary  conditions  of  the 
wing  surface  are  taken  directly  from  wing  surface  pressure 
distributions  measured  by  experiment,  that  is, 

Vm  2'  (4) 


We  select  the  exterior  boundaries  of  pressure  measurement  control 
surfaces  close  to  the  top  and  bottom  walls  of  the  wind  tunnels.  We 
still  use  the  experimentally  measured  pressure  distributions  to  be  the 
boundary  conditions.  If,  during  wind  tunnel  experiments,  static 
pressure  holes  are  opened  to  measure  pressure  on  the  side  walls  of  the 
upper  and  lower  flows  of  the  test  section,  then,  the  boundary 
conditions  for  the  upper  and  lower  flows  can  be  directly  taken  as 
experimental  pressure  distributions.  If  wind  tunnel  structure  is  a 
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limiting  factor,  and  there  is  no  way  to  make  experimental  measurements 
of  the  pressure  distribution  on  the  boundaries  of  the  upper  and  lower 
flows,  then,  the  boundary  conditions  for  the  upper  and  lower  flows  of 
the  wind  tunnels  can  be  approximated  by  making  use  of  the  already 
known  conditions  at  the  four  corner  points  of  the  area  of  calculation 
to  combine  them  into  a  multi-term  expression  and  obtain  the  required 
values  f9l*  In  this  way,  the  boundary  conditions  for  all  the 
calculations  of  flow  fields  inside  wind  tunnels  are  Dirichlet  boundary 
conditions. 


2.  Boundary  Conditions  for  Calculations  of  Free  Flow  Fields 


As  is  shown  in  Fig.  2,  wing  surface  boundary  conditions  for  the 
calculation  of  free  flow  flow  fields  are  taken  from  experimental 
pressure  distributions  or  already  given  pressure  distributions. 


<**<>#»«  =  # 


u-0£iif»0  j 


Fig.  2  Boundary  Conditions  from  Calculations  of  Free  Flow  Flow  Fields 
1 .  Or  2.  Airfoil 


As  far  as  the  boundary  conditions  of  distant  fields  are  concerned, 
from  the  basic  physical  nature  of  the  problem,  it  is  possible  to  know 
that  at  an  infinitely  distant  point  m«o  •  If,  during 
calculations,  the  distant  field  boundaries  are  selected  sufficiently 
far  from  the  model,  it  is  possible  to  take  the  boundary  conditions  for 
an  infinitely  far  point  and  make  an  adequate  approximation  of  the 
boundary  line  [101.  During  the  calculations  in  this  article,  we 
selected  the  upper  and  lower  flow  distant  field  boundaries  at  a 
distance  10-15  times  the  chord  length  of  the  model.  V/hen  considering 
the  special  point  from  which  transonic  horizontal  disturbances 
propagate,  the  upper  and  lower  areas  should  be  even  farther  from  the 
airfoil  at  approximately  20-25  times  the  chord  length. 
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IV.  Mixed  Difference  Forms 

This  article  makes  use  of  non-equidistant  difference  forms.  The 
second  order  derivative  in  the  y  direction  makes  use  of  the  central 
difference 


2  ht/i-i  U..H- (Ay,.,+  Ay,)  M.-.rf  Ay, 
Ay,-,  Ay,  (Ay,.,  +  Ay,) 


(5) 


The  derivative  in  the  *  direction  must  depend  on  the  velocity 
discriminant  for  that  point 


1  —  =  -V+ 1.  Mi  u 


to  determine  its  form. 

When  (i  —\I:) >o  ,  "-  and  u '««  select  for  use  the  central 

difference 


u 


si.  1“ 


j.J  “M.  -Jj. 

Ax,.,  + Ax, 


i.,=* 2 


A x, .,  «, », j -  <  A x,.,  +  Ax,)m,.,4-Ax,  m,-i. l 
Ax,.,  Ax,(Ax,.,  + Ax,) 


(7) 

(8) 


When  u’-A/'XO  ,  «.  and  “«»  adopts  for  use  the 

difference  on  one  side  of  the  upper  flow 


u. 


i.i 


u , —  u,  -i.i 

Ax,..  +  Ax,., 


u„  ,  ,  -  2  (Ax,.,+  Ax.-Jm,.,,^  Ax,.,  u, i 

Ax,.,  Ax1.t(Ax,.l  +  Ax,-l) 


(9) 

I 


(10) 


When  ,u -  A/-')  -o,  ,  the  equation  degenerates.  However,  the  », 

term  still  adopts  for  use  the  central  difference  form  (7). 

This  article,  in  its  treatment  of  shock  wave  points,  selects  for 
use  a  "shock  wave  capturing  method."  According  to  a  special 
characteristic  of  the  pressure  equation  velocity  discriminant  (6), 
when  (i-A/1), and  u  -i/'),.j>o  f  a  shock  wave  step 
develops  between  two  points.  Choosing  the  points  ft  as  shock 
wave  points,  reference.  [9]  carried  out  numerical  value  tests  of  three 
types  of  forms.  It  was  discovered  that  taking  shock  wave  points  and 
comparing  them  to  entering  subsonic  velocity  points  is  appropriate. 

V.  The  Linear  Relaxation  Iterative  Substitution  Solution  Method  /1 62 

In  turn,  at  the  point  (i  ,  »  according  to  the  form  of  the 

equation  for  that  point,  we  select  differing  difference  forms  and 
substitute  in  equation  (2).  At  the  same  time,  we  take  the  boundary 
conditions  and  insert  them  into  the  difference  equations.  Because  of 
this,  it  is  possible  to  obtain  a  set  of  high  order  non-linear 
algebraic  equations  relating  to  »  .  This  article  uses  a  linear 

relaxation  iterative  substitution  method  of  solution,  choosing  the 
line  of  relaxation  along  the  v  direction. 

According  to  the  numerical  value  tests  carried  out  in  reference 
£ 9 j ,  if  it  is  desired  to  make  the  calculations  for  subsonic  and 
supersonic  flow  movements  all  converge,  then,  the  term  u\  must  be 
used  in  the  calculation  of  the  earlier  field,  that  is, 


\  A  / 

( <V’ v 

V  A'v,_,  -!  \r,_.  ) 


(1  -A/1  ),.;>() 
a  -  ihf  |< o 


(ID 


If  we  use  the  newest  values  in  the  calculations,  it  will  lead  to 
divergence  in  the  calculation  of  supercritical  flow  movements. 
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This  article  makes  use  of  the  previous  field  value  to  calculate 
velocity  discriminants  16J.  The  parameter  for  the  term  adopts 

the  use  of  earlier  field  values  in  calculations.  The  remainder  all 
make  use  of  the  newest  values  in  calculations.  After  rearrangement, 
one  finally  obtains  a  one  by  three  diagonal  matrix 

Oi  !  'V',  1,-rf,  (12) 

Equation  (12)  uses  a  "pursuit"  type  of  method  for  its  solution.  Each 
time  calculations  are  done,  they  begin  by  choosing  a  zero  initial  flow 
field.  In  order  to  speed  up  convergence  as  well  as  improving  the 
stability  of  convergence,  this  article  makes  use  of  relaxation 
calculations.  The  selection  of  relaxation  factors,  on  the  basis  of 
the  calculation  experience  in  reference  [9],  can  make  use  of 
superrelaxation  in  the  case  of  purely  subsonic  critical  flow  movements. 

As  far  as  supercritical  flow  movements  are  concerned,  the  use  of 
superrelaxation  makes  it  easy  to  cause  divergence  in  calculations. 

The  adoption  of  a  low  relaxation  is  relatively  more  appropriate.  WhenA/„S;0.9 
one  should  adopt  super  low  relaxation  factors  (for 
example,  -c»n*0.3)  ).  Only  then  is  one  able  to  guarantee 

stable  convergence  in  numerical  value  solutions. 

VI.  Calculation  Tests 

In  order  to  test  the  feasibility  of  the  methods  discussed  above, 
this  article  used  a  velocity  potential  method  sample  calculation  to 
carry  out  comparative  calculations  with  reference  [5]  and  reference 
[6].  Table  3  shows  the  calculated  results  of  pressure  distribution  in 
this  article  and  in  reference  [5]  for  the  point  on  an  airfoil  i.u 
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Fig.  3  Pressure  Distribution  for  the  Point  >’/c»o.u  on  the  Airfoil 
'M. » o.78 1. <i*2*56')  l.  Calculated  Results  With  This  Article's  Metod 
2.  Calculated  Results  from  Raference  [5]_ 


/1 63 

The  airfoil  is  a  supercritical  airfoil,  v./u.i ,  at'  •  In 

the  Pig. ,  the  results  for  this  article  are  obtained  by  solving  small 
perturbation  velocity  potential  equations  on  the  basis  of  interior  and 
exterior  boundary  conditions  from  test  pressure  distributions. 


Pig.  4  Calculated  Free  Plow  Pressure  Distribution  At  Point  •*,c“  ‘*PS 

i.  Results  Calculated  i.n  Referen.ee  [6]  2.  Results 

Calculated  in  This  Article 
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Pig*  5  Calculated  Free  Plow  Pressure  Distributions  At  Point 
1 .  Calculated  Results  from  Reference  [6]  2  Calculated  Results  from 

This  Article  aa,  *o.9,««o*> 

Pig.  4  and  Pig-  5  present  the  curves  comparing  the  calculated 
results  at  the  point  y/c- l.os  for  the  pressure  distribution  of  free 
flow  flow  fields,  obtained  with  the  wing  surface  pressure  distribution 
as  inner  and  outer  boundary  conditions  using  the  NACA  0012  airfoil  in 
reference  [6],a-0°,  A/„  =  o.72and  0.9  *  Prom  Pig.  5-5  it  is  possible  to 

see  that  the  results  calculated  in  this  article  and  in  [5]  and  [6] 
agree  relatively  well.  This,  therefore,  demonstrates  that  the 
convergence  solution  achieved  with  the  method  in  this  article  is 
reliable. 

In  order  to  check  out  a  step  further  the  convergence  speed  of  the 
calculation  method  in  this  article,  this  article  ran  a  comparison  of 
the  convergence  speeds  of  the  mixed  difference  method  for  solving 
pressure  equations  [9]  and  for  solving  velocity  potential  equations 
for  large  disturbances  in  the  x  direction  [6].  Due  to  the  fact  the 
factors  actually  influencing  iterative  substitution  convergence  speed 
are  very  numerous-including  such  factors  as  the  grid,  the  initial 
field,  the  accuracy  of  convergence,  and  the  actual  boundary 
conditions-it  is  only  possible  to  carry  out  a  rough  comparison  of  the 
two.  The  results  of  this  are  shown  in  Table  1  . 
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Table  1  Comparison  of  Convergence  Speeds  in  Numerical  Value  Solutions 
for  Pressure  Distribution  Equations  and  Velocity  Potential  Equations 
(The  Airfoil  Used  Is  NACA  0012, a=0*  )  2.  Item  3*  Method  4- 

Solution  of  Pressure  Equation  [9]  5*  Solution  of  Velocity  Potential 

Equation  [6]  for  Large  Disturbances  in  the  '*  Direction  6.  Grid  7. 
Number  of  Iterative  Substitutions  8.  Time  9.  Siemens 


Prom  the  Table  it  is  possible  to  see  that,  when  the  method  in  /164 
this  article  and  the  method  in  reference  [6]  are  compared,  the  method 
in  this  article  shows  relatively  faster  iterative  substitution 
convergence  and  uses  relatively  less  time. 


VII.  Typical  Sample  Calculations  for  Aspects  of  Applications 

When  the  author  of  this  article  did  research  on  the  problems  of 
tunnel  wall  interference  in  transonic  airfoil  wind  tunnels,  he 
produced  the  concept  of  using  a  mixed  difference  method  to  directly 
solve  equation  (2).  After  awaiting  numerical  value  tests  of  its 
stability  and  convergence  characteristics,  he  then  turned  the  tables 
and  used  it  for  calculations  of  tunnel  wall  interference  and  in  the 
calculation  of  difficult  problems  involving  solving  for  airfoil 
configurations  for  already  known  pressure  distributions.  Its  basic 
principles  and  details  can  be  seen  in  reference  [9]*  Typical  sample 
calculations  can  be  seen  in  Pig.  6-9. 
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Fig.  6  Evaluation  of  Calculated  P . c  vr.vv  distributions  on  Control 
Surfaces  <y/e«i.oa>  naca  ooi2,«*o\<»*=2?s  1.  Free  Flow 

Evaluation,  Calculated  Values  2.  Test  Values 

Fig.  6  presents  a  comparison  of  the  calculated  pressure 
distribution  and  the  empirical  pressure  distribution  in  an  evaluation 
of  free  flow  in  this  article.  This  comparison  is  made  for  a  NACA  0012 
airfoil  in  a  situation  where  M  =0.6,  0.75,  0.9>  a- 0%  and  the 
open-closed  ratio  cr««  2  %  .  Because  of  this,  it  is  possible  to  make 

a  definite  determination  of  the  size  and  nature  of  the  amount  of 
tunnel  wall  interference.  In  the  Fig.,  the  distance  from  the  control 
surface  to  the  model  is  u/c=  1.05  ». 

Fig.  7-8  respectively  show  the  pressure  distribution  when  we  use 
the  circular  arc  given  in  [ll]  for  M  =0.707  aad  0.817,  , 

followed  by  the  results  of  a  solving  back  for  the  geometrical 
configuration  of  the  airfoil. 

Fig.  9  shows  the  airfoil  configurations  designed  from  the 
pressure  distributions  which  came  out  of  this  article  for  the  NACA  0012 
airfoil  presented  in  [6],  when  M  ■  o.85,  a==o°  •  From  Fig. 

7-9  one  can  see  that  the  airfoil  configurations  calculated  in  this 
ar.ticle  match  up  quite  well  when  compared  with  precisely  determined 
airfoil  configurations. 
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Pig.  7(a)  Experimentally  Determined  Pressure  Distribution  on  Airfoil 
Surface  (11)  >M-»u.ro7.  «  =  <>•> 


Fig.  7(b)  Airfoil  Configuration  Calculated  from  Empirical  Pressure 
Disribution  1 .  Already  Known  Airfoil  2.  Airfoil  Obtained  from 
Calculations  of  the  Experimentally  Measured  Pressure  Distributions  in 
7(a)  3  Circular  Arc  Airfoil 


Pig.  8(a)  Experimentally  Determined  Pressure  Distribution  on  Airfoil 
Surface  (11)  *M~»o.*i7,  «-<>•> 
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Pig.  8(b)  Airfoil  Configuration  Calculated  from  Experimentally 
Determined  Pressures  1  .  Already  Known  Airfoil  2.  Airfoil  Obtained  by 
Calculations  from  the  Experimentally  Determined  Pressure  Distribution 
in  8(a)  3.  Circular  Arc  Airfoil 


VIII.  Conclusion 

The  preliminary  research  which  this  article  has  done  on  numerical 
value  calculation  methods  for  transonic  small  disturbance  pressure 
equations  has  already  demonstrated  that,  on  condition  that  an 
appropriate  difference  form  is  selected,  as  well  as  an  appropriate 
iterative  substitution  method  and  relaxation  parameter,  the  use  of  a 
mixed  difference  method  to  solve  equation  (2)  produces  convergence  and 
that  its  convergence  solutions  are  reliable. 

While  satisfying  the  assumed  condition  of  small  disturbances, 
this  method  handles  tunnel  wall  interference  problems  with  pressure 
distribution  boundary  conditions.  The  method  in  this  article  handles 
boundary  conditions  with  ease,  converges  quickly,  and  has  such 
advantages  as  saving  machine  time.  Moreover,  the  solutions  which  are 
derived  for  flow  field  pressure  distributions  can  also  directly  act  as 
the  basis  for  comparing  definite  tunnel  wall  disturbances,  avoiding 
transformations,  being  both  simple  and  direct.  Besides  this,  the 
method  in  this  article  is  used  in  designing  airfoil  configurations 


# 


from  given  pressure  distributions.  It  has  also  obtained  relatively 
satisfactory  preliminary  results. 

Research  into  the  calculation  methods  and  applications  of  small 


disturbance  pressure  equations  has  still  just  begun.  There  are  still 
a  good  many  questions  awaiting  in-depth  research  and  solution. 
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Rig.  9(a)  NACA  0012  Airfoil  Pressure  Distribution  (6)  Mft.ss,  a*or 


Fig*  9(b)  Airfoil  Configurations  Calculated  from  Given  Pressure 
Distribution  1.  Exact  Airfoil  2.  Airfoil  Obtained  from  Calculations 
of  Experimentally  Determined  Pressure  Distributions  9(a) 
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NUMERICAL  SOLUTION  OF  TRANSONIC  SMALL 
DISTURBANCE  PRESSURE  EQUATION  USING 
A  MIXED  DIFFERENCE  NETHOD 

Wang  Lixia  Luo  Shijun 
(Northwestern  Polytechnical  University) 

Abstract 

A  transonic  small  disturbance  pressure  equation  (i.  e.  TSDP  equation) 

u)u.,+  «M-£±lA/d«i-0 

is  proposed  for  computing  transonic  flow  fields  in  wind  tunnel  or  free 
streams. 

The  mixed  difference  method  is  used  to  calculating  TSDP  equation. 
Numerical  experimentation  indicates  that  the  use  of  suitublc  difference 
schemes  and  relaxation  technique  yields  converged  solutions  to  TSDP  equa¬ 
tion.  Comparisons  show  that  TSDP  solutions  agree  well  with  those  of 
transonic  small  disturbance  potential  equation.  Application  of  the  proce¬ 
dure  to  assessing  transonic  wind  tunnel  interference  and  designing  airfoil 
from  the  givcu  pressure  distribution  arc  illustrated. 
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SUMMARY  Using  the  point  of  view  of  vortical  dynamics,  we  have 
taken  a  new  look  at  the  mutual  effects  of  moving  bodies  and  fluids, 
and  we  have  obtained  an  even  deeper  understanding  of  this  classic, 
fundamental  problem.  This  article  does  research  on  the  effects  of 
bodies  on  vorticity  fields-  and  gives  a  generalized  theory  of  the 
production  of  vorticity  on  the  surface  of  objects  as  well  as  fo.r  its 
dissipation  in  fluids.  It  also  presents  two  types  of  vorticity 
sources  on  the  surface-  of  objects-the  total  and  the  local-and  analyses 
the  effects  of  each  of  the  sources 

I.  Introduction 

One  of  the  fundamental  problems  of  fluid  dynamics  is  the  mutual 
effects  of  moving  bodies  and  fluids.  This  includes  both  the  aspect  of 
the  effects  of  bodies  on  fluids  and  the  aspect  of  the 
counter-reactions  of  fluids  to  bodies.  The  traditional  route  for 
handling  this  problem  is  through  a  Navier-Stokes  equation  analysis  of 
the  mutual  relationships  between  the  velocity  field  v  and  forces 
received  by  the  surface  of  the  objects.  However,  this  process  alone 
is  not  able  to  exhaust  the  knowledge  relating  to  the  mutual  effects  of 
fluids  and  moving  bodies.  In  order  to  check  out  a  new  path,  we  looked 
at  the  considerations  below. 

The  special  point  with  fluids  is  that  stress  tensors  are  not 
directly  related  to  either  velocity  fields  or  velocity  gradient  fields  -vf  . 

As  is  universally  known,  Vu  can  be  uniquely  analysed 
into  a  vorticity  field  ,  an  expansion  field 

and  a  irrotational  non-diffusion  field.  For  real  fluids,  vorticity 
fields  are  restrained  by  the  object  surface  viscosity  conditions,  and 
completely  irrotational  non-diffusion  fields  are  almost  non-existent. 

(1  ]  tA)  Because  of  this,  the  process  of  moving  bodies  affecting 
fluids  is  actually  the  process  of  objects  producing  vorticity  fields 
and  expansion  fields.  Conversely,  the  situation  of  forces  being 
received  by  objects  is  also  not  the  same  velocity  field  itself.  It  is 
very  closely  related  to  the  two  derived  fields  of  vorticity  and 
expansion.  Because  of  this,  it  is  possible  to  predict  that  directly 
researching  the  mutual  effects  of  these  two  derived  fields  and  objects 
will  enrich  and  deepen  the  knowledge  we  obtain  through  the  use  of  the 
traditional  method  forming  a  second  generation  of  theory  on  the  mutual 
effects  of  fluids  and  bodies. 


As  far  as  non-compressible  fluids  are  concerned,  there  is  only  a 
single  vorticity  field.  The  unique  mechanism  of  bodies  transfer ing 
energy  to  non-compressible  fluids  is  none  other  than  the  production  of 
vorticity  through  adhesion  boundaries.  Because  of  this,  the  reaction 
of  fluids  to  bodies  is  also  none  other  than  the  effects  of  the  already 
produced  vorticity  on  the  bodies.  The  object  of  this  research  is,  on 
the  basis  of  this  type  of  simple  typical  situation,  to  set  up  a- 
gener&l  theory  of  the  mutual  interactions  of  vorticity  fields  and 
moving  bodies. 

This  research  is  composed  of  two  parts.  These  respectively 
consider  the  two  sides  of  the  mutual  effects  discussed  above.  Part  I 


(the  main  text)  sets  up  a  general  theory  for  the  production  of 
vorticity  from  any  given  movement  over  any  given  surface  as  well  as 
for  its  dissipation  in  the  fluid.  The  leader  in  this  kind  of  research 
is  Lighthill  (3]>  Recently,  because  of  the  fact  that  research  /169 

concerning  separation  flows  and  vortical-  movements  has  been  extensivly 
developed,  he  has  also  come  in  for  renewed  attention  (Reference  l4]  in 
several  reports  from  the  1935  AIAA  Conference  on  Shear  Plow  Control.) 

This  article  develops  the  ideas  of  Lighthill  into  a  finished  form 
giving  a  scalar  expression  for  vorticity  sources.  This  enables  us  to 
make  a  comprehensive  and  systematic  solution  for  the  mechanism 
producing  vorticity.  These  results  have  actually  provided  several  new 
theoretical  foundations  for  the  processes  controlling  shear  layers  on 
the  surface  of  objects  and  their  separation. 


II.  A  General  Formula  for  the  Strength  of  Vorticity  Sources 


When  the  forces  penetrating  the  bodies  are  strong,  we  take  the 
vortical  dynamics  equation  for  compressible  fluids  and  integrate  the 
fluid  volume  v,  •  Moreover,  we  set  up,  except  for  the  object  surface 
on  the  fluid  boundary.  It  is  possible  to  demonstrate 

that 


— -j  £  n  -V5  dV 
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II 


is  normal  to  the  outside  of  the  surface  of  the  object. 
Because  of  this,  ^  is  a  unique  vorticity  source.  Moreover, 

o  =*  —  »  n  •  V  ta 

is  the  amount  of  vorticity  produced  in  a  unit  time  for  a  unit  surface; 
that  is,  the  vorticity  source  strength  as  first  introduced  by 
Lighthill  [3] .  In  the  case  of  a  two  dimensional  flat  plate  with 
movement  of  uniform  speed,  through  a  simple  proof,  reference  [31 
actually  gives  us: 


—jr 


on 


dB 


(2) 


Now,  we  generalize  these  results  to  the  effects  of  any  given  movement 
on  any  given  curved  surface  and  derive  a  general  formula  pulling 
together  the  mutual  relationships  between  3  >  p  and  the  surface 

friction  *  .  We  often  need  normal  gradients  and  tangential 

gradients.  In  order  to  provide  these,  it  is  best  to  introduce  the 
dynamic  coordinate  framework  (5),  which  is  often  used  on  the  curv 
in  differential  geometry.  The  vector  radius  for  any  point  on  S  fron 
the  fixed  origin  point  can  be  written  as  £•<«*  =  1,2)  and 

is  a  parameter.  Moreover, 


a; 

i)s* 


(  r  ,  x  r  i#0) 


that  is,  the  two  basic  stress  vectors  forsthe  point  in  question  on  the 
tangential  plane  n  It  is  also  possible,  from  to 

introduce  the  basic  contravariant  vector  ^  *  as  well  *  -  the  measuring 
tensor  0.s=r.*r,  or  g^-r'-r1  .  Let  the  determinant  of  9 . 

be  the  unit  normal  vector  of  u,  ,  the  unit  normal  '^sctor  n  of  •£ 
is  then 


and 


Therefore,  {?,}  —  Cf ,,  rr,  J)  and  {r*}  =  (rl,  r1,  n) 

respectively  form  the  three  dimensional  stress  and  contravariant  right 


hand  dynamic  coordinate  frame 


<  r  =1,  2,  J), 


on 


■S’  as  is  shown  in 


Pig.  1  .  It  is  obvious  that  the  determinant  of  the  corresponding  three 
dimensional  measurement  tensor  j=1.2,3)  is  still  g  . 


If  we  set  up  as  the  variable  along 


then,  it  is  obvious  that 


the  normal  derivatives  of  r.  and  n  are  r,.,« 0,  n\,--o. 

Their  tangential  derivatives  from  Gauss  formulas  and  Weingarten 
formulas  are 


|  r..a  /"A  ri+b.tn,  r‘. - /"V*  r 

n.i®-6S  r.~.r-6.,r“  (3) 

In  these  results,  r.\=»r and  is  the  two  dimensional  second 

type  Christofell  symbol  on  5  .  -  7  .  7$  is  the  s 

second  basic  equation  parameter.  Here  and  below,  we  see  repeated 
examples  of  the  customary  practices  of  differential  geometry  such  as 
taking  basic  vectors  and  writing  them  into  calculation  formulas.  We 
will  frequently  take  a  certain  vector  /  and  operation  V  and 

resolve  them  into  vector  projections  in  the  normal  direction  on  the 
tangent  plane  *  .  At  this  time,  we  write 

7~/*  i  /  •  ~/*  r  .»  7>~n/,=  n/’ 

V(*)  —  V.(*)i  n  (•).,«=  r  ’  (•),.+  b  (•). , 


Besides  this,  let 
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This  is  a  component  of  a  third  order  totally  skew  matrix.  Here, 
contains  the  familiar  three  index  symbols.  We  then  have 


r"‘  ■ c , 


n  x  /  -if’*'  f.r  , 


Through  the  use  of  this  equation  and  (3),  it  is  possible  to  verify  the 
identity 

«  x  <  /  *Vn  -  /V-  «)=»-(«  >  /  >  •  V  »  (4) 


It  is  now  possible  to  derive  a  general  expression  for  <j  . 
First  of  all,  solve  for  the  normal  gradient  of  «  on  any  given  curbed 
surface  5  in  the  fluid.  After  that,  change  it  to  a  form  suitable 
for  use  on  the  surface  of  the  object  dB  .  Set  the  viscous  stress  on 
to  be  7  We  have 

7  =#i(A  u  +Vt/>  •  n  (5) 

is  the  transpositon  of  Ay  ,  and  it  is  easy  to  prove  that 


7  =*<( a  x  n  +  a  > 


In  this 


^  AM  ^  *• 

A  ”2V  y  •  n 


See  reference  [6]  (According  to  reference  [6]  a  is  2?v.«,  •  There 
seems  to  be  an  error.)  After  that,  we  note 


We  have 


/i&3.=  n  *  ff  **  n  x  (7  ~t*A)  (3) 
Because  of  the  fact  that  n.i  =  0,  f  we  obtain 

n  x§  ,i+p n  a *, 

In  this,  on  the  basis  of  V-5=o  ,  we  have 

©*,  -4—  (g) 

*  ^ 


Therefore,  from  (8),  we  obtain 

Here,  we  make  use  of  V» x  «  = ^  *  r  '  =o » 

Due  to  the  fact  that  the  normal  derivative  of  'B  breaks  off 
during  the  period  when  it  penetrates  dB  ,  it  is  necessary  to  make 
some  additional  changes  to  "  x5.i  .  From  (8)  we  have 

^VXu.-fi-Vn-D.)+"V-B-BV-n 

Therefore,  by  the  use  of  identity  (4)»  it  is  possible  to  obtain 

n  XB:> "  »  x  (Vx  a)  +i*n  x(VXb,)-(S  X  jj)  ?V  n 

In  this,  it  is  possible  to  show  that 

n  X  (V  x  ®  »>  *V.  ®. 

Besides  this,  the  Navier-Stokes  equation!  gives  us 


-mx(vx5) 


;x(v„+p!l») 
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By  combining  the  results  above,  we  obtain  no  for  any  given  curve  on  S’ 


M».i-  n  x  (vp  +  P^)-(n  x  g  >  -V  n  +  n  C  n  •  (9  X  fl  )]  +  m  V.  o>,  ^  ^ 

Here,  we  have  already  avoided  any  interruptions  in  normal  gradients  on  OB 
How,  we  take  S.to  be  the  object  surface  OB,  .  On  it,  the 
velocity  at  any  point  is 


;»i/(i)+fl(()x'r'  (11) 

o  ia  the  translation.  £  is  angular  velocity.  These  are  only 
functions  of  t  .  ?'*?-?.  is  the  vector  radius  which 

corresponds  to  the  instantaneous  center  of  rotation.  r»  is  the  vector 
radius  for  that  instantaneous  center  in  a.  fixed  coordinate  system. 
Because  of  the  fact  that  ,  r,.s=r..0  .  On  the  basis 

of  adhesion  conditions,  ( 1 1 )  is  nothing  other  than  the  fluid  point  dB 
mass  velocity.  If  we  take  ( 1 1 )  and  substitute  it  in  (6),  we  obtain 

'  Cl)  •  (  r  ,  x  n  )  j  on  dB 


Therefore, 


A  A. 


,  and,  it  is  possible  to  obtain 


Thi3  also  gives 


Therefore 


n*A=-2Q.  (12,.) 

A  «  »  X  (a  X  n  )  *2  n  X  & ,  (l  -  j) 


V%A  n  —  fl.  V*  n  )>  n,(Vx^)_0 


In  another  respect,  we  have 


V.«,=2  V.(fl*n)  =  -26.,  ; 


—  2i5 « *  V  n 


In  this  way,  if  we  take  (7)  and  substitute  it  into  the  rigot  side  of 
(10),  it  is  possible  to  discover  that  the  term  A '  along  ao  n  (iV.ar 
disappear  together.  Therefore,  if  we  cake  7  on  cue  surface  of  the 
object  and  write  it  as  r »  ,  we  obtain  a  general  expression  for  ° 


3=4 


(-nXVp  +  (nxf).v?-»7Cn-(VX?)]}-3x^- 

a/ 


(13) 


In  this,  dv  Idt  is  easy  to  obtain  from  (11). 

Particularly  in  the  case  of  two  dimensional  flows,  we  have 

i-W+t!)  (h) 


Reference  [4]  has  already  used  (14)  to  discuss  the  mechanism  for  the 
production  of  vorticity  in  two  dimensional  smooth  flow  on  flat  plates 
for  different  acceleration  movements.  Reference  [7]  makes  use  of  (14) 
and  a  simplified  point  vortex  system  model  to  make  a  generalized 
discussion  concerning  non-steady  state  separation  in  two  dimensional 
uniform  flows.  It  also  makes  definitive  solutions  for  a  series  of 
experimental  results  concerning  non-steady  state  shear  layer  control. 

In  another  regard,  when  ,  due  to  'one  fact  that  r-*o,  , 

even  in  the  case  of  three  dimensional  flow,  we  still  have  sgu -.-.cion  /172 

(14).  Lighthill  [3]  has  already  pointed  out  that,  at  this  time,  it  is 
not  n-*o,t  hut,  it  is  '3. .  The  vorticity  is  concentrated  i  AB 
the  unbounded  thin  layer  and  is  not  dispersed  toward  the  interior  of 
the  fluid.  Here,  one  should  remember  that  the  right  side  of  { 1 ^ ) 
comes  from  the  viscosity  term  of  the  NS  equation  and  adhesion 
conditions.  Reference  [4]  only  determines,  on  the  basis  of  the  surface 
form  for  (14),  that  the  source  of  3  is  non-viscous.  This  thesis  is 
not  sound. 


The  derivation  of  (13)  is  used  in  movement  equations  ana 
viscosity  condltons.  It  is  also  unobtrusively  used  in  the  continuity 
equation  v?=0,’  •  Due  to  the  fact  that  the  district  t ;  ;  .  o  pt  ? 

on  dB  is  not  known,  (12)  -is  not  the  formula  for  caleuU. i(J  o  .  oh; 
is  the  expression  which  relates  the  distributions  of  a  .*;u  />,  f 

well  as  accelerations  on  the  surface  of  the  object.  Attention  should 
be  paid  to  the  assumption  that  nt  can  be  differentiated.  Because  of 
this,  it  is  required  that  the  surface  of  the  object  be  everywhere 
smooth. 

III.  The  Mechanism  Producing  Vorticity  on  the  Surface  of  Objects 

Prom  (13)  one  can  clearly  see  that,  on  the  surface  of  objects, 
there  are  two  types  of  vorticity  sources-? 

(I)  Stress  Sources 

These  types  of  sources  exist  for  any  movement  on  any  curved 
surface.  For  certain  types  of  classical  body  forms,  they  have  special 
characteristics  which  can  also  be  called  local  sources.  These  are 
al30  composed  of  three  parts. 

1.  Tangential  Source  5„=-(l/p)  nxv/> 

This  type  of  source  has  already  been  discussed  in  references  L3], 
[4],  and  [7].  Just  as  is  pointed  out  in  [7],  if  it  is  possible  to 
ignore  the  pressure  gradient  which  is  caused  by  boundary  layer 
displacraent  thickness  on  flat  plates  by  semi-limitless  smooth  flows  of 
uniform  velocity,  then,  it  is  only  in  the  vicinity  of  the  forward  edge 
stationary  point  that  o,,*0, 

2.  Tangential  Source  Vi,*(1/p)(b  x  t  )  •  V  if 

This  is  an  effect  of  object  surface  curvature  and  the  three 
dimensional  coupling  oi  Iron  ( 0 )  and  (12b),  it  is  also  possible 

to  write  this  as 


(15) 


Because  of  this,  except  for  rotation  of  the  body,  at  the  point 
the  w  which  the  object  surface  already  h-.?  ..Ltl  bring  about 
the  uninterrupted  x^roduction  of  new  vortic.u.y  in  fluids.  On  sharp 
edges  where  Vn  is  extremely  large..  °»»  is  very  important.  At  such 
points,  o.p  ie  also  very  large.  These  two  together  determine  the 
form  of  free  vortical  layers  on  sharp  edges.  Discussion  below  points 
out  that  a  *•  determines  the  effects  of  friction  forces  on  the 
object  surface  as  a  whole. 

3.  Normal  Source  a,,  -  n  (r>~  -  (  n  /p>C  n  *(V  x  f  )] 

Prom  (9)  ,  it  is  m’.so  x->ossible  to  write  this  as 

5.,=w.-(3.  (16) 

Because  of  this,'0,1  -px-^ars  in  rotational  areas  of  the  two  dimensional 
vector  field:?  jn  the  object  surface  or  in  areas  of  dissipation  in  the  , 
field.  For  example,  various  areas  on  the  rotating  object.  In 
the  case  of  non-rotating  objects,  the  discussion  below  applies. 

(II)  Acceleration  Sources 

When  objects  make  non-inert ial  movements,  on  the  one  hand,  it 
will  influence  the  />»  ?  distribution  and  alter  the  stress  source. 

(15)  presents  an  example  of  this.  On  the  other  hand,  it  also  produces 
the  tangential  source  -tixduldt,  which  depends  only  on  movements 
as  affected  by  surface  geometry.  This  is  recorded  as  a .  Idea 
(11)  we  have 

5  *#--=  “  «  x{f)  +fjx  r  '  +  _Q  x  fj'  f  x  (Q  x  r  ) } 

(17) 


In  two  dimensions 


CT, 


X  0  -/}(»•?')  -Q  (  £  .Q  )  {-(„* 


)D’ 


(IS) 


The  simplest  example  is  simple  translation  acceleration  movements  [4]. 
Consider  another  type  of  example  of  the  two  dimensional  oscillation 
of  wing  edges  as  shown  in  Fig.  2.  Using  (13)  and  the  symbols  in  the 
Fig.  it  is  easy  to  demonstrate  that,  on  the  top  and  bottom  surfaces, 


one  has 
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o  ,m*=e  ,r  Q\e  i-ti  x  e  i,  r=  j  f 


Note  that  this  effect  has  no  relationship  to  the  direction  of 
rotation.  Moreover,  on  the  trailing  edge..  (AB  surface),  one  then  has 


°'*  =  e  ;  CQ 


is  the  chord  length. 


Pig.  2  Two  Dimensional  Oscillating  Wing  Edges  or  Flaps 


It  should  he  pointed  out  that  surface  vorticity  a  and  vorticity 
sources  are  two  different  concepts.  According  to  Taylor's  development 
of  this,  o  only  applies  to  vorticity  in  the  interior  of  designated 
fluids.  This  is  in  agreement  with  equation  (1).  On  the  other  hand, 
from  (8)  and  (12b)  we  have 


5  =  a  (0ndB) 


(13a) 


f  )  *  n 


(19b) 


Fig.  3  Velocity  Forms,  Vorticity  Forms,  and  Vorticity  Sources  (a) 
Velocity  Forms  (b)  Vorticity  Forms  (c)  Vector  Vorticity  Source 
Distribution  (d)  Scalar  Vorticity  Source  Distribution  / 1 74 

Because  of  this,  object  surface  u  is  directly  determined  by  f  .ud  ; 

Moreover,  it  is  not  produced  by  <5  .  Particularly,  in  the  case 

of  the  two  dimensional  encircling  flows  when  SJ'0  ,  at  times  when 
the  coordinate  system  is  fixed  to  the  surface  of  objects,  one  has 


Moreover 


o  —  vn  X  v)  )  —  —  v  n  x  5  . 


Because  of  this,  <3  and  3  are  respectively  as  shown  by  the 
velocity  forms  and  the  rates  of  slope  and  curvature  on  the  surface  of 
objects  in  Pig.  3* 

If  a  body,  beginning  at  instant  ,  begins  to  move  in  a  fluid 

which  was  originally  stationary,  stress  sources  ana  acceleration 

sources  will  both  begin  to  come  into  play.  in  particular, 

during  the  process  of  initial  movement,  will  start  to  exert  an  effect 

over  the  whole  surface  of  the  object.  When  the  object  trausi fcions  to 

movements  of  uniform  speed  at  time  /<  ,  one  has  for  t*'t> 

*  • 

It  is  only  the  remaining  local  sources  °  and  .r  , 

which  continue,  on  portions  of  the  object  surface,  to  send 
vorticity  toward  the  interior  of  the  flu-id.  The  amount  which  is  sent 
is  closely  related  to  the  <3  distribution  on  the  surface  involved 
Because  of  this,  global  acceleration  sources  and  localized  stress 
sources  have  very  different  natures. 


IV.  Dissipation  of  Scalar  Sources  and  Vorticity  in  Fluids 


Fig.  3  shows  that,  when  «  and  6  on  the  surface  of  the  object 
are  opposed  to  each  other,  the  velocity  form  becomes  concave.  This 
foretells  the  coming  of  flow  separation.  Because  of  this,  the 
relationship  between  the  directions  of  -  and  6  are  an  important 
basis  for  the  determination  of  separation.  In  order  to  do  this,  it  is 
necessary  to  consider  the  scalar  source  o  •  a  -  -  a/2)  m  •Vo* 

This  quantity  makes  it  possible  for  us  to  distinguish  sources  and 
confluences.  If  we  consider  the  situation  in  which  the  acceleration 
form  is  convex,  (u»-o>0>  is  caxied  a  source.  Then,  exists  a 

situation  in  which  w  •  a  <0  is  a  confluence.  Therefore,,  the 
appearance  of  confluences  and  separation  are  closely  related. 


Taking  the  situation  in  which  -M  "  0  as  an  example,  if  we  make 
use  of  (13)*  (17),  and  (19b),  we  obtain 


at 


a  —  — 


1  - 

>  y 


1  V p  +  vat  •  V  n  •  at  —  v 


(20) 


The  contribution  from  normal  sources  disappears.  If  we  let  r >C,  . 
then, ' obviously ,  the  back  pressure  gradient  and  the  acceleration 
movement  both  produce  vorticity  confluences.  In  particular,  from  (3). 
we  have 


it  *V  n  •  a  1 —  fc . »  o’  *  (a* 

The  second  basic  rearrangement  of  this  on  a  curved  surface  L 5 ( number 
unclear)]  is 


v.  b^dfdV 

To  look  at  it  another  way,  because  of  the  fact  that  the  symbols  in 
are  only  related  to  the  geometrical  form  of  the  surface  of  the 
object,  this  has  potential  usefulness  for  controling  flow  separation. 

o<t)  -  (  *  at1  dr 

*  1 

In  the  same  way  that  the  vector  source  *  functions  in  (l), 
scalar  sources  influence  the  characteristics  of  integration 
operations.  This  integration  is  called  enstrophy.  The  scalar 
quantity  V  here  should  not  be  confused  with  the  angular  velocity  of 
solid  body  rotation  vector  sj  previously  introduced.  If  we  integrate 
after  taking  the  dot  product  of  o>  and  vorticity  dynamics  equations, 
we  obtain  [8] 


dQ 

dt 


Lj  ( w  *D«  w  -j.  V  w  «Voj)  d\' 


VJ  2~  "  *  Vw'  dS 


(21  ) 
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In  this,  D  is  the  rate  of  strain  tensor.  *  represents  a  second 
order  contraction  of  the  tensor.  The  last  term  in  equation  (21)  is 
the  overall  scalar  source.  In  equation  (21),  stands  for 

the  enstrophy  density  transformation  created  because  of  pulls  and 
bending  forces  on  vortex  tubes.  When  two  dimensional,  this  is  zero. 
Putting  this  in  an  orthogonal  curved  coordinate  system 
causes,  for  any  (x,D.x,  to  be  always  along  the 

direction.  Then,  « *D* Z  -cj* D,,  .  Because  of  this,  the 

existence  of  the  mean  value  /),,</)  35(1/2)  Day  causes  us  to 

obtain 

.  j  fi)*D •<adl'  -  DiOilU) 

*t 


In  another  regard, 


<Pmv 


J 


V<3  (U  dh'^o 


This  is  the  attenuation  of  vorticity  fields  caused  due  to  viscosity 
dissipation.  Because  of  this,  we  write  the  overall  scalar  source  as  If. 
and  solve  (21)  for  the  enstrophy  energy  operational  formula 


QU) 


(22) 


For  two  dimensional  flow,  we  simplify  and  get 


(0)  +j  '  [^(r) -0(r)]Jr  (23) 

It  is  obvious  that,  in  (22),  the  functions  of  index  factors  and  of 
are  different.  Only  the  latter  reflects  the  creation  and 
disappearance  of  OU)  .  There  is  a  popular  point  of  view  which 
recognizes  that  viscosity  will  not  cause  vorticity  to  disappear  but 
will  simply  cause  it  to  be  redistributed.  This  is  incorrect.  In 
actuality,  if,  at  a  certain  instant,  the  moving  body  is  taken  out, 
then,  it  is  plain  to  see  from  (23)  that  &  will  ultimately  attenuate 


to  zero.  Therefore, 


6  ss  o 


V.  Conclusions 

(I)  Vorticity  in  noncompressible  fluids  is  produced  by  object 
surfaces  through  adhesion  conditions.  It  is  diffused  by  viscosity  and 
dissipated  by  it. 

(II)  On  the  surface  of  moving  objects,  there  are  two  different 
types  of  vorticity  sources.  One  type  is  an  overall  or  global 
acceleration  source  related  to  movements  and  object  surface  geometry. 
The  other  type  is  the  localized  stress  source.  This  includes 
tangential  sources  created  by  pressure  gradients,  tangential  sources 
from  object  surface  curvature  and  vorticity  coupling,  as  well  as 
normal  sources  of  rotation  areas  in  r  vector  fields.  The  latter  two 
types  of  localized  sources  only  appear  in  three  dimensional  flows. 
These  localized  source  areas  are  areas  of  energy  supply  from  which 
energy  is  transfered  to  fluids  by  objects  in  uniform  motion 

(ill)  Except  for  vorticity  and  vector  sources,  it  is  possible  to 
introduce  scalar  sources  as  the  effects  of  area  sources  and 
confluences.  The  existence  of  vorticity  confluences  is  the  first  sign 
of  flow  separation.  The  effect  of  three  dimensional  curved  object 
surfaces  on  scalar  sources  is  actually  determined  only  by  the 
characteristics  of  object  surface  geometry  and  has  potential 
usefulness  in  the  control  of  flow  separation. 
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INCOMPRESSIBLE  THEORY  OF  THE  INTERACTION  BETWEEN 
MOVING  BODIES  AND  VORTICITY  FIELD 

- THE  GENERATION  OF  VORTICITY  BY  BODY 

SURFACES  AND  ITS  DISSIPATION 

Wu  Jiezhi 

(Chinese  Aeronautical  Establishment) 

Abstract 

The  interaction  between  moving  bodies  and  fluids,  n  classical  problem 
of  fluid  dynamics,  is  reexamined  from  the  viewpoint  of  vorticity  dynamics. 
In  this  way,  we  may  gain  some  new  insight  into  the  mechanism  of  the 
interaction  and  can  be  led  to  a  series  of  results  which  arc  of  practical 
value.  The  present  paper  studies  the  action  of  a  moving  surface  to  vorti¬ 
city  field  and  gives  a  general  incompressible  theory  of  the  generation  of 
vorticity  at  the  surface  and  its  dissipation  in  the  fluid.  It  is  found  that 
there  are  two  types  of  vorticity  sources,  the  global  one  depends  only  on 
the  acceleration  property  of  the  surface  geometry,  while  the  local  one 
exists  in  both  accelerated  and  uniform  motion,  consisting  of  tangential 
sources  from  pressure  gradient  and  a  three-dimensional  effect  of  the  sur¬ 
face  curvature,  and  a  normal  source  due  to  the  divergence  of  the  two- 
dimensional  vorticity  on  the  surface. 


A  STRONG  INV I SC ID- VISCOUS  INTERACTION  SOLUTION 
OP  A  PLANE  TRANSONIC  CASCADE  PLOW 
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Chen  Yunwen  Shen  Mengyu 
(Qinghua  University) 


Zhang  Yaoke 

(Computing  Center,  Academica  Sinica) 


I.  Introduction 


In  order  to  design  transonic  turbine  wheel  machinery  for  both 
high  .efficiency  and  high  loads,  in  aerodynamic  calculations,  it  is 
necessary  to  deal  with  the  calculation  of  the  influence  of  viscosity. 
When  gas  movements  on  the  surface  of  turbine  blades  do  not  produce  se 
paration,  or  when  areas  of  separation  are  very  small,  gas  movements  in 
the  flow  way  can  be  divided  into  turbine  blade  surface  and  trailing 
edge  lower  flow  thin  viscous  layer-separation  layers  and  wake  areas  as 
well  as  non-viscous  flow  areas  other  than  these.  Because  of  this, 
it  is  possible  to  distinguish  the  flow  movements  in  these  two  types  of 
areas  which  should  be  solved  for  by  the  use  of  Euler  equation  sets  and 
those  which  should  be  solved  for  with  boundary  layer  equation  sets. 

The  mutual  influences  of  iuviscid  and  viscous  flows  become  apparent 
through  a  system  of  iterative  substitution  calculations  of  inviscid 
flows  and  boundary  layers.  By  carrying  out  inviscid  calculations  from 
given  aerodynamic  parameters  and  blade  cascade  geometric 
configurations,  we  obtain  a  flow  movement  parameter  distribution  for 
turbine  blades  and  wake  surfaces.  This  acts  as  input  data  for  t tie 
initial  calculation  of  boundary  layers  and  wakes  to  obtain  the 
exclusion  thickness  •  for  boundary  layers  and  wakes.  Moreover,  we 
take  this  geometric  boundary  of  inviscid  flow  and  carry  out 
adjustments  so  that  it  becomes  the  flow  boundary  for  the  next  inviscid 
calculations .  This  type  of  iterative  substitution  is  repeatedly 
carried  out,  and  it  continues  straight  on  until  the  complete  body 
convergence  conditions  for  the  whole  field  are  satisfied,  and  then  it 
stops.  The  flow  chart  for  these  iterative  substitutions  is  as  shown 
in  ?ig.  1  . 

In  these  calculations,  we  make  use  of  a  time  advance  integration 
method  for  finite  areas  to  calculate  inviscid  flow  fields  111,  We 
make  use  of  an  integration  relationship  equation  to  solve  for  laminar 
flow  and  turbulent  flow  boundary  layers  as  well  as  wakes  We 

make  use  of  empirical  relationships  to  precisely  determine  the  turning 
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point  position  ( 4 ] •  After  all  the  calculations  converge,  on  the  basis 
of  [5]»  we  carry  out  calculations  of  the  mixing  of  gas  flows  after  the 
cascade.  Prom  this,  we  obtain  the  loss  coefficient  for  the  blade 
cascade.  Below,  we  give  a  more  detailed  description  of  the  carrying, 
out  of  the  matters  mentioned  above.  Finally,  we  give  two  sample 
calculations  for  the  flow  movements  in  transonic  turbine  wheel  flat 
plane  blade  cascades. 


Pig.  1  1  .  Determining  of  the  geometrical  configuration  of  the  blade 

cascade  and  other  aerodynamic  and  geometrical  parameters  for  intake 

and  exhaust  ports.  2.  Inviscid  flow  calculations  5*  Is  there 

convergence  or  not?  4«  Boundary  layer  and  wake  calculations  5* 

Alterations  of  flow  field  boundaries  6.  Yes  7  Plow  mixing 
co.lculation3.  Calculations  of  blade  cascade  loss  coefficient  8. 

Machine  stop 


II.  Inviscid  Plow  Calculations 

If  we  assume  constant  specific  heat  in  a  perfect  gas  during 
adiabatic  flow,  we  make  use  of  a  time  advance  method  to  solve  for  its 
steady  state  movements,  and,  at  this  time,  the  basic  equation  set  for 
the  planar  movements  can  be  taken  to  be: 
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(2.2) 


dS*  +  ~£t  <*>+/>“">  +  ~-<P»v)  =  0 

*dF  +  -^T(f)uv)+  (p  +  pv’)^0  (2.3) 

*  -f  +  T(w,  +  y,)=W  (2.4) 

Here,  />  are,  respectively,  the  density  and  pressure  of  the  g-.  . 
are  the  velocity  components  of  the  gas  along  the  y  directions-  c> 
is  the  specific  heat  at  fixed  pressure.  R  is  one  gas  constant. 
is  the  overall  heat  content  or  enthalpy  of  the  gas  at  the 
intake.  As  is  shown  in  the  Fig.,  ABCDEFGH  is  the  area  of 
calculation.  GPF  and  BSC  respectively  the  equivalent  boundaries  of  t> 
pressure  surface  and  suction  surface  of  the  two  phases  sho/  i  <ir-  ' 
blades.  The  intake  boundary  and  the  exhaust  boundary  are 

parallel  to  the  v>-  axis.  The  parallelogram  ABGH  is  the  periodic 
flow  field  area  in  front  of  blades.  The  included  angle  between  the 
straight  lines  ,AB.  and  HG  and  x  axis  are  the  intake  flow  angl: 

When  making  initial  inviscid  cal  cud  ti.ons,  FE,  Cr> ;  are 
the  directions  along  the  estimated  exhaust  angle  0,  .  It  later 
inviscid  calculations,  CD  ^re  basically  wake  boundary  lines 

solved  for  in  initial  wake  calculations,  and,  in  the  process  of 
calculations ,  continuously  undergo  adjustments. 

In  the  ;T  direction,  it  is  possible  to  take  equal  distances  (or 
unequal  distances)  and  make  contour  lines  roughly  parallel  to  the  y 
axis.  In  the  flow  path,  it  is  possible  to  take  equal  distances  (or 
unequal  distances)  and  make  lines  giving  a  rough  simulation  of  the 
flow.  The  points  at  which  tne  contour  lines  or  interval  lines 
intersect  with  the  flow  simulation  lines  are,  then,  the  grid  points. 

The  points  of  actual  calculations  are  located  between  mutually 
adjacent  grid  points.  This  is  represented  in  Fig.  3  with  the  symbol.  • 


We  take  the  grid  unit  of  two  points  between  which  there  is  a 
centrally  located  calculation  point  and  form  from  it  a  unit  of 
calculation.  If  we  take  the  equations  (2.1 )— (2-3)  after  they  have 
been  made  non-dimensional,  and  change  the  unit  for  these  calculations 
to  an  integral  form,  it  is  possible  to  obtain: 


|  j  dy  +  (£  (pa,  pu)«n</s  =  0 

li  dx  dy  +  (p  + p  u‘,  puv)»nds  —  0 

ii  d.x  dy  +  <j)  (puv,  p  +  pv1)*  n  ds<~ 0 


(2.6) 


(2.7) 


Here,  .  n  is  an  exterior  normal  line  unit'vector  for  the  boundary  of  a 
calculation  unit.  According  to  Denton  [6],  the  concept  of  setting 


cr-  ■—  m  i.  ■ .  i  a 
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mm 
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up  a  difference  grid  is  that,  carrying  out  a  dispersion  of 
(2.5)— (2.7) *  it  is  possible  to  obtain  the  apparent  calculation  form. 
From  the  flow  parameters  for  the  time  increment  n  ,  it  is  possible  to 
calculate  the  flow  movement  parameters  for  the  time  increment  n+l 
The  form  of  its  calculation  expresses  the  equations  refereu  go  in  ilj 
and  [ 7 ] . 
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The  calculation  boundary  conditions  are  as  shown  below.  The 
intake  and  exhaust  boundaries  AH  An.l  DE,  when  axial  velocities  are  0<«<c 
(  c  is  the  speed  of  sound),  on  .AH  ,  we  are  given  u.V\. 

as  well  as  the  gas  entry  angle  ift,  ■,  Cn  DE  ,  w?  re  given  the  back 
pressure  p,»:  on  blade  surface  use  M-s  *'•?’-!  the  mutual  intersection  oC 
the  equivalent  object  surface  boundaries  and  inviscid  flow  speeds  on  CPF 
On  All  and  HG  ?  ,v,3  select  for  use  periodicity  conditions  for 
flow  move  :r ■ v  :  to  act  as  conditions  for  determining  solutions. 

CD,  fE 

are  wake  boundary  lines.  "Based  on  the  exclusion  functioning  of 
wakes  in  inviscid  flows  (main  flows)  as  well  as  the  equivalence  of  gas 
static  pressures  along  wake  cross  sections,  it  is,  therefore,  true 
during  calculations,  that,  after  each  time  increment,  at  corresponding 
points  on  CD  and  Ft  ,  the  pressure  takes  its  average  value,  in  order 
to  guarantee  pressure  equivalence  at  corresponding  points.  At  the 
same  time,  on  CD  and  FE  ,  we  eliminate  the  normal  velocity  component 
of  gas  flow  in  order  to  guarantee  contact  between  flow  movement 
velocities  and  wake  boundaries.  The  positions  of  CD  n  FE 
determined  by  moving  the  wake  exclusion  thickness  on  either  side  ox 
the  wake  base  line.  When  doing  preliminary  boundary  layer 
calculations,  we  take  the  wake  base  line  to  be  the  flow  line  passing 
through  the  wake  edge  points.  In  the  later  iterative  subsitution 
calculations,  the  curve  passing  through  wake  edge  points  and 
contacting  the  corresponding  points  on  the  two  sides  of  the  wake  in 
the  average  direction  of  the  gas  flow  is  taken  as  the  wake  base  line. 

III.  Boundary  Layer  and  Wake  Calculations 

1  .  Laminar  flow  boundary  layer  calculatio  ns. 

We  first  select  for  use  the  Illingsworth-Stevartson 
transformation.  We  take  the  differential  equation  set  for 
compressible  laminar  flow  boundary  layers  on  adiabatic  object  surfaces 
and  transform  them  into  the  differential  equation  sets  for 
non-compress ible  laminar  flow  boundary  layers.  After  this,  we  make 
use  of  the  single  parameter  Lot tsianski i  method  to  solve  for  its 
momentum  integration  relationships.  Finally,  through 
re transformation ,  one  solves  for  and  obtains  the  kinetic  energy  loss 


thickness,  the  exclusion  thickness,  as  well  as  the  wall  surface 
coefficient  of  friction  (Refer  to  reference  [2]  )  for  compressible 
laminar  flow  boundary  layers. 


2.  Determination  of  Turning  Points 

According  to  the  nature  of  turning  points,  it  is  possible  to 
divide  them  into  natural  turns  and  gas  bubble  turns.  In  these 
calculations,  if  gas  bubble  turns  are  produced,  let  us  assume  that  the 
gas  bubble  length  is  zero.  That  is,  let  us  recognize  that  laminar 
flow  separation  points  are  the  starting  points  for  turbulent  flow 
boundary  layers.  If  natural  turnings  are  produced,  we  need  make  no 
assumptions  about  the  length  of  the  turning  area.  In  the  second 
situation,  the  boundary  layer  kinetic  energy  loss  thickness  after 
adopting  the  turning  point  does  not  change.  The  critical  Reynolds 
number  and  boundary  layer  thickness  ratios  after  turning  points  are 
gotten  from  empirical  formulas  [4]. 

3.  Turbulent  Plow  Boundary  Layer  and  Turbulent  Plow  Wake 
Calculations 

In  these  calculations,  we  select  for  use  the  Lag-Entrainment 
integration  relationship  method  of  such  people  as  Green  [3]  to 
calculate  turbulent  flow  boundary  layers  and  turbulent  flow  wakes. 

The  Lag-Entrainment  method  of  such  people  as  Green,  through  the 
generalizations  of  such  people  as  East  [8],  is  used  in  the  calculation 
of  turbulent  flow  boundary  layers  in  areas  having  small  separations. 

At  certain  places  before  the  separation  points,  the  calculations  turn 
toward  an  "inverse  form",  avoiding  the  appearance  of  a  singularity  in 
the  original  equation  set  at  the  separation  point.  At  this  time,  the 
boundary  layer  exclusion  thickness  «M.r<  is  a  given  value. 

Moreover,  the  inviscid  flow  velocity  »r(x)  ls  a  value  awaiting 
determination.  In  the  calculations,  it  is  ^.eessary  to  make  (180) 

continuous  adjust  n?nts  to  the  given  d.(.v),  /  causing  the  calculated 

values  of  end  obtained  from  inviscid  flow 

■ -lculr  'on  ;■  to  coincide. 
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In  the  calculations,  it  appears  that,  in  the  case  of  certain 
transonic  turbine  blade  cascades  (such  as  sample  calculation  1),  the 
gases,  while  passing  through  the  blades,  receive  an  abrupt 
acceleration,  leading,  in  the  vicinity  of  the  rear  edge  of  the  blades, 
to  a  change  in  flow  parameters  exceeding  the  appropriate  ranges  of  bhe 
radial  empirical  formulas  given  by  such  people  as  Green.  This  makes 
it  impossible  to  carry  out  the  calculations.  At  this  time,  use  is 
made  of  the  IT  as  h-McD  o  nal  d  method  [9]  which  is  introduced  inLo  the 
calculations  after  the  reactions  3top  in  order  to  carry  them  out. 

This  overcomes  the  difficulties  mentioned  above. 

IV.  Iterative  Substitution  Calculations  and  Convergence  Principles 

1.  First  of  all,  we  carry  out  time  advance  inviscid  calculations 
for  blade  pathways  of  a  given  geometrical  configuration.  This  leads 
directly  to  the  convergence  of  the  calculations,  and  the  convergence 
principles  are: 

<i)  E,  =  Max( | <«;,*;  - i /p  > <0.0002 

H 

or 

-^-KIC'-Cil/FXO. 00005 

Here 

k-[^-  E  <«■+»■>::: ]'" 

N  is  the  total  number  of  calculation  points.  '-^present  the  rank 

and  file  designations  of  the  calculation  points.  The  superscripts  !l 
and.  ”  '  1  are  the  number  of  time  increments. 

2.  From  the  results  of  inviscid  calculations,  one  calculates  the 
exclusion  thicknesses  of  boundary  layers  and  wakes  Moreover,  one 
carries  out  corrections  of  blades  and  wake  boundaries.  In  order  to 
improve  the  stability  of  calculations,  one  carries  out,  on  boundary 
layer  thicknesses,  hyporelaxed  iterative  substitution: 


0<6)<1 


0  is  the  relaxation  parameter, 

3.  The  results  of  the  invito l--.  Calculatiions  above  act  as  an 
initial  calculation  field.  If  one  carries  out  a  certain  incremental 
inviscid  calculation  on  the  adjusted  flow  path,  and  one  arrives  at  the 
overall  iterative  substitution  convergence  principles  below,  then  the 
calculations  stop.  If  not,  then,  the  iterative  substitution  continues 
to  progress. 

The  overall  iterative  substitution  convergence  principles  are: 

(1)  The  calculation  errors  between  consecutive  time  increments  in 
inviscid  flow  calculations  are  the  same .- (4 . 1  ) 

(2)  The  maximum  error  between  two  adjacent  ,Af  number 
distributions  on  blade  surfaces  in  inviscid  calculations  (also  the 
input  data  for  two  adjacent  boundary  layer  calculations)  is 

Max  |  M  "’-A/,"’!  <0.005 

1 

(3)  The  maximum  difference  value  between  the  blade  surface 

coordinate  •  !/,.,.{  in  inviscid  boundary  calculations  and  the 

coordinates  for  the  equivalent  blade  forms  obtained  in  the  same 
boundary  layer  calculation  is 

Max |y, •.<«!-!/<. nl ;  cascade  distance  ,<> 

1 

The  selection  of  the  value  of  &l  should  be  related  to  the 
precision  which  it  is  possible  to  achieve  in  the  •i-ich!  ning  of  blades. 
In  the  sample  calculations,  we  select  4  =0.0003,  which  far  exceeds 
the  actual  blade  machining  precision. 


1 


V.  Calculations  of  Losses  and  Gras  Exit  Angles 
After  overall  convergence  of  iterative  substitution,  according  to 
the  Stewart  theory  [5],  other  aerodynamic  and  geometrical  parameters 
for  the  blade  cascade  as  well  as  the  boundary  layer  exclusion 
thickness  for  blade  rear  edge  suction  surfaces  and  pressure  surfaces 
and  their  associated  kinetic  energy  loss  thickness  are  used  to 
calculate  blade  cascade  loss  parameters  and  mixing  effects  of  blade 
cascade  lower  flows  as-  they  impact  on  wakes  and  main  flows.  At  the 
same  time,  on  the  ba3is  of  the  direction  of  tne  wake  base  line  after 
the  convergence  of  the  cal  .culations,  it  is  possible  to  make  precise 
determinations  of  numerical  values  for  gas  exit  angles. 


VI.  Sample  Calculations 


Example  one.  This  is  a  calculation  of  the  middle  section  of  a 
blade  cascade  in  a  624  turbine.  The  principle  aerodynamic  parameters 
are  as  follows.  Gas  entry  angle  ^,  =  50.4%  .  Overall  entry 


A  =“295987. 3  N/M- 


Overall  entry  temperature  7\«=  288  K 


presure 

Back  pressure  P*~  102730.3  N (M1  „ 

Example  two.  Inis  is  an  RA  turbine  blade.  Principle  aerodynamic 
parameters  are  as  follows.  Gas  entry  angle  £.=60%  .  Overall 

entry  pressure  />«  =  230300.3  N/M1,  .  Overall  entry  temperature 

71  =  288  K,  .  Back  pressure  Pi ~  101292.8  NJ M‘  » 

The  number  of  grid  points  for  example  caluclation  number  one  and 
example  calculation  number  two  were  respectively  41x11  and  71x11. 

After  convergence  of  the  inviscid  calculations,  we  carried  out 
boundary  layer  calculations  and  adjusted  the  inviscid  calculation 
boundaries.  After  that,  after  inviscid  calculations  were  carried  out 
following  each  50  time  increments,  one  boundary  layer  calculation  was 
carried  out,  and  this  continued  right  on  until  we  arrived  at  overall 
convergence.  Fig.  4  is  the  example  calculation  one  blade  surface 
non-dimensional  pressure  coefficient  -P-PlPa>  -calculation 
results  and  experimental  values  (  °*  are  respectively  the 

critical  density  and  the  critical  sonic  speed  corresponding  to  entry 
gas  flow  parameters.)  Pig.  5  i°  the  number  distribution  and 
experimental  results  for  the  sample  OHl.culaT.ion  two  blade  surfaces. 

In  this  Pig.,  we  also  see  presented  the  calculated  results  for 
inviscid  flow. 
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Pig.  4  1.  Inviacid  2.  Viscous  Flow-Inviscid  Iterative  Sues ti tat  ion 

3*  Experimental  Results 

Pig.  5  1 .  Inviscid  2.  Viscous  Flow-Inviscid  Iterative  Substitution 

3.  Experimental  Results 


The  degree  of  correspondence  between  inviscid-viscous  flow 
iterative  substitution  calculated  results  and  experimental  values  i3 
an  obvious  improvement  over  inviscid  flow  calculated  results.  Besides 
this,  we  also  calculated  the  kinetic  energy  loss  para-lets'* .  In  the 
case  of  the  624  central  blade  cascade,  its  v-Aue  s  £  **  0.04921, 

.  Prom  the  wake  base  line  direction,  we  determined  the  gas  exit  angle 
to  be  -6(5.75%  •  These  values  are  in  relatively  ju 

agreement  witn  the  corresponding  experimental  values  /?,=  -G7.4° 

(taken  from  [lO]).  In  the  case  of  the  RA  turbine  blad.- 

cascade,  the  calculated  example  for  £  is  0.0248.  .  Because  there 

are  no  experimental  results,  there  is  10  vay  to  make  a  comparison. 


VII.  Several  Observations 


1 .  Calculations  demonstrate  that  the  positions  of  wake  boun- 


daries  and  the  specifying  conditions  for  the  wake  boundaries  have  a 
very  great  effect  on  calculated  results  for  blade  surface  aerodynamic 
parameters.  We  went  through  the  process  of  carrying  out  test 
calculations  on  the- 624  middle  blade  cascade  for  the  two  respective  ■ 
situations  below-  When  the  wake  is  fixed  on  the  flow  lines  passing 
through  the  edge  of  the  wake  as  obtained  in  initial  inviscid 
calculations,  as  well  as  for  the  case  where  periodic  boundary 
conditions  are  taken  as  the  defining  conditions  for  wake  boundaries, 
the  differences  between  blade  surface  pressure  distributions  and 
inviscid  flow  calculation  results  are  extremely  small. 

2.  Trailing  edge  thickness  of  blades  has  a  very  large  influence 
on  blade  cascade  loss  coefficients.  This  is  particularly  true  of 
blades  which  have  relatively  large  trailing  edge  thicknesses.  The 
error  in  loss  parameters  caused  by  errors  in  blade  thicknesses  will  be 
nuch  larger  than  the  errors  caused  by  calculation  errors  from  boundary 
layer  thicknesses.  Because  of  this,  when  one  carries  out  comparisons 
of  calculated  values  and  experimental  results,  it  is  necessary  to 
carry  out  caluclations  on  the  basis  of  the  actual  thicknesses  of  the 
experimental  blade  models  used. 

3.  What  this  article  describes  is  the  initial  results  of  our 
work.  Looking  from  the  standpoint  of  the  calculation  results  in 
the  preliminary  work  we  have  just  done,  and,  using  the 
inviscid-viscous  flow  iterative  substitution  method  presented  in 
this  article,  the  degree  of  correspondence  between  the  results 
obtained  and  experimental  values  is  quite  good.  However,  this 
work  still  waits  for  even  more  examples  of  calculations, 
including  the  carrying  out  of  more  sophisticated  calculations  and 
checks  of  flow  movements  on  blades  having  small  separation  areas. 
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A  STRONG  INVISCID-VISCOUS  INTERACTION  SOLUTION 
OF  A  PLANE  TRANSONIC  CASCADE  FLOW 

Chen  Yunwen  Zhang  Yaoke  Shen  Mengyu 

0 Qinghua  University)  ( Computing  Center,  Academia  Sinica)  (Qinghua  University) 

Abstract 

An  inviscid-viscous  interaction  method  has  been  developed  to  predict- 
the  bludc-to-blade  flow  in  plane  turbine  cascades.  The  interaction  effect 
is  taken  into  account  by  iteratively  solving  the  inviscid  and  viscous  flo.ws. 
The  inviscid  flow  is  calculated  by  a  finite  urea  time-marching  method. 
The  viscous  flow  is  calculated  by  the  integral  method,  Loitsianskii  s  meth¬ 
od  for  laminar  boundary  layer  after  performing  Illingworth-Stewartson 
transformation  and  Green's  lag-entrainment  method  for  turbulent  boundary- 
layer  and  wake.  A  mixing  calculation  is  carried  out  to  determine  the 
kinetic  energy  loss  coefficient.  The  exit  angle  is  determined  by  the  base 
line  of  the  wake  predictions  agree  well  with  the  test  results  for  two 
transonic  turbine  cascades. 


NUMERICAL  CALCULATIONS  OF  INVISCID  TRANSONIC 
FLOWS  OVER  WINGS 
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Chen  Zuobin  Yao  Furu  Zhang  Yulun 


(China  Aerodynamic  Research  and  Development  Center) 


SUMMARY  This  article  introduces  a  type  of  numerical  value  method 
of  calculation  for  non-viscous,  steady  state  transonic  flows  around 
wings.  It  selects  for  use  the  full  potential  equation  as  mathematical 
model.  Through  the-  use  of  appropriate  coordinate  transf ormations , 
sweptback  wings  are  changed  to  rectangular  ones.  Also,  an  unlimited 
physical  domain  is  changed  to  a  limited  calculation  domain.  In  the 
domain  of  calculations,  use  is  made  of  the  mixture-type  limited' 
difference  form  to  discretize  equations.  The  set  of  algebraic 
equations  formed  from  difference  equations  is  solved  through  linear 
relaxation  iterative  substitution.  Successive  refinements  of  the 
calculation  grid  make  the  calculations  relatively  economical. 


I.  Introduction 

Wings  are  the  aircraft's  most  important  aerodynaiuo  force 
components.  The  quality  of  their  aerodynamic  characteristics  has  an 
extremely  great  influence  on  the  capabilities  of  the.  whole  aircraft. 
Moreover,  the  development  of  research  into  the  calculation  of 
transonic  flows  around  wings  not  only  has  significance  for  academic 
knowledge.  It  also  has  very  great  applied  value. 

The  numerical  value  method  presented  in  this  article  is  the 
solution  for  transonic  speed  full  potential  equations  in  an  orthogonal 
coordinate  system.  We  first  use  shear  transformations  to  change 
backswept  wings  into  rectangular  wings.  After  that,  we  use 
compression  transformations  to  change  an  unlimited  physical  domain 
into  a  limited  domain  of  calculation.  In  this  way,  on  the  one.  hand, 
we  can  make  direct  use  of  distant  field  boundary  conditions  and,  on 
the  other  hand,  it  is  also  possible  to  get  an  automatic  concentration 
of  grid  points  in  the  vicinity  of  the  wings  and  make  a  reasonable 
distribution  of  the  grid  points.  The  coordinates  after  the 
transformations  are  still  orthogonal.  The  treatment  of  boundaries  and 
relaxation  scanning  are  directly  observed.  This  article  is  a 
continuation  and  a  development  of  the  work  in  reference  £  1 ] . 
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When  the  surfaces  of  objects  satisfy  precise  boundary  conditions, 
the  actual  method  for  handling  them  is  similar  to  the  wing  surface 
boundary  treatment  in  Carlson  airfoil  calculations  [2].  In  every 
cross  section  of  wings,  we  take  the  perturbation  speed  potential 
analysis  and  expand  it  to  the  imaginary  grid  point  simulation  inside 
the  wings  in  order  to  use  it  to  set  up  difference  equations  for  grid 
points  close  to  the  surface  of  objects-  At  points  of  subsonic  speed, 
we  select  for  use  central  difference  forms.  At  points  of  supersonic 
speed,  we  select  for  use  rotational  difference  forms  as  presented  by 
Jameson  [3].  This  is  done  in  order  to  guarantee  that  the  difference 
forms  will  have  accurate  areas  of  dependence. 


II.  Basic  Equations  and  Boundary  Conditions 


Under  conditions  where  there  are  no  strong  shock  waves  and  steady 
state,  non-viscous  compressible  gas  movements,  it  is  possible  to  make 
an  irrotational  assumption,  and  satisfy  the  full  potential  equation: 

(aJ  —  u!)<px,+  (a1  —  v')<p„  -b  (a:  —  w:)<p!a  —  2uv  q>,,-2vvj  — 2 >nv  <P„  =  0  < j  ) 

In  order  to  eliminate  the  singularity  at  an  infinitely  far  point,  we 
introduce  the  perturbation  speed  potential 

G  —  <f-  x  cos  a  —  y  siti  a  (2) 

From  velocity  potential  equation  (1),  it  is  possible  to  deduce  that 
the  perturbation  speed  potential  G  satisfies  the  equation  below: 

1  >35 

(u*  u*)G tt  +  (a1  - u'K?,,  +  (a1  - u>')(7„ - 2uv  G„ - 2 uw  G„ -  2vw  G,,~ 0  f  3  '1 

Velocity  components  are 
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u  =  Videos  u-i-G, 
v  =95, —  sin  a  +  G , 
w=Vj  =  G, 


(4 


The  local  speed  of  sound  is 


r 

The  pressure  parameter  #  is 

c'"5Bri[1+  (i> 

All  the  physical  quantities  in  equations  (l)-(6)  are 
non-dirnensionalized  with  the  basic  chord  cu  an(^  the  incoming  flow 
velocity  g 

Object  surface  boundary  conditons:  velocity  and  object  surfaces 
make  mutual  contact.  Say  that  the  upper  and  lower  wing  surface 
equations  are 


y 


!/*(*.*> 


upper  surface 


lower  surface 


(7) 


Then,  the  object  surface  boundary  conditions  can  be  represented  as 

u?9et 


,  ,  uPPe* 

(cos  a  4  G,)  — —  (sina+-(7,)  +  G,-^-  —  0 


dz 


(cos  a 


+C.,^-La+C.,+G.^< 


(6) 
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Vortical  surface  conditions:  ignore  the  initiation  of  curling  of  wake 
vortical  surfaces.  Also,  assume  the  wake  vortical  surface  is 
congruent  with  the  3(02  coordinate  surface.  Then,  on  the  wase 
vortical  surface 

I  'Jj 

m 


Distant  field  boundary  conditions: 


C“0  (jc=  -oo,  y,z=*  ±co)  v  1 

£»»+Cr„=0  (*  =  -f  oo) 

(12 


Symmetry  conditions:  on  the  plane  of  symmetry  2-0  conditions 
should  satisfy 

W-G'= 0 


III.  Coordinate  Transf ormations 

In  order  to  make  calculations  convenient,  we  make  a  rational 
arrangement  of  the  grid  points,  first  transforming  sweptback  wings 
into  rectangular  wings.  The  actual  transformation  is  as  follows: 

X  «•  <x  -  xLg  (z) )  /cU)  -  0. 5 

K-tf 
Z~  z 


After  that,  we  take  a  limitless  physical  domain  and  transform  it  into 
a  limited  domain  of  calculation.  In  this  way,  one  is  then  able  to 
directly  use  distant  field  boundary  conditions.  One  first  selects  an 
appropriate  transformation  function,  which  is  still  able  to  make  most 
of  the  grid  points  fall  close  to  the  surface  of  the  wing.  The  actual 
transformation  is  as  shown  below: 

In  the  x  direction  and  the  y  direction 
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1 36 


*  -< 


-x.  +  At  ig  |  'l  (:  |  •!-  A,  tg  j  ^  {..)  |  X<,  -  X. 

-Xr~X::Xt 

•v,  '  At  tr  f  l  ■  X.\  j  :  .i,  ,c  j  *  C  -  j  A- A', 
U=At  tg(  7  j 


(15) 


(lo) 


The  ‘  direction  transformation  is  similar  to  the  direction 


transformation.  The  transformation  is  seen  in  Pig.  1 . 


Fig.  1(a)  The  Relationship  Between  the  Physical  Coordinate-  Cyc tea  -x.y.z 
and  the  Calculation  Coordinate  System  1  .  Pl*ne 

Fig.  1(h)  The  Relationship  Between  the  Physical  Coordinate  Sye  •  r  \  x.v.z 
and  the  Calculation  Coordinate  System 
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After  the  coordinate  transformations,  equation  (3)  becomes 


f C''(*Gz),+DfgG.'+EfhGi.+ 

+  F9ftG,:+Gf(fl/Gol  +  Hfh(/&\):=o  *  (17) 


The  velocity  components  (4)  become 

u  ~co$  a  +  /GK/c(z) 

'  v  =*sin  a-\  gG,  ,'3) 

w  ~  / >  I  <r  i+h  G t 


IV.  Difference  Forms 

At  points  of  subsonic  speed,  central  difference  forms  are  always 
selected  for  use  At  points  of  supersonic  speed,  in  order  to  produce 
a  reasonable,  man-made  viscosity,  and,  to  make  the  difference  forms 
stable,  we  make  use  of  localized  velocity  coordinate  systems.  Then, 

( l~  ■fr)<P»+(Av,-v’..)  =  o  (^3) 


In  the  calculation  coordinate  system 

f{f  GOi+Bt  g(g  G\),  +C,  h(h  C?c);  +2D,  fgGt,+ 
+  £,  fh  GiL+Fl  gh  G,:+Gt  /(/.  /  C<) ,  +  /?,  /A(/,  C?t) .] 

Av5=  c(V)  +  /AG,.+ 

+  /./</, /G’t),+//l(/lCt); 


92 


In  order  to  cause  the  difference  equations  to  have  precise  areas  of 
dependence,  at  points  of  supersonic  speed,  we  select  for  use  a 
rotational  difference  formula  [3j*  That  is,  in  the  term  {*$>-*>„)  f 
the  derivatives  of  the  various  terms  use  a  central  difference.  In  the 
term  .  the  various  term  derivatives  use  the  headwind  difference. 

The  difference  equation  set  is  solved  "by  the  use  of  linear  relaxation 
iterative  substitution.  The  iterative  substitution  process  can  be 
seen  a's  being  the  introduction  of  an  artificial  time  coordinate  term  A/, 
Because  of  this,  in  the  equations,  except  for  *•••  'f>mt  , 

we  also  added  the  terms  **«»»  $>„,  .  -Jameson  makes  use  of  the 

aggregate  difference  between  the  initial  iterative  substitution  value 
and  the  last  iterative  substitution  value.  Let  W»»  he 

implied  in  the  difference  formulae.  Jameson  has  already  pointed  out 
that,  in  order  to  guarantee  precise  stability,  it  is  sometimes 
necessary  to  add  the  term  V>u  .  This  method  is,  in  a  real  form, 
adding  v,i»  <p„,  implied  in  -.  Because  of  the  fact  that, 

when  there  is  convergence,  v„ .  <r.<  tends  toward  zero* 

therefore,  in  this  way,  there  is  no  influence  on  the  final  results. 

The  actual  difference  forms,  when  r>0,w>o  ,  are  exemplified 

by  the  various  terms  in  <p.. 

(sG%)n  ™ .*)  —  g,.  l(G{,l.l.t—G,.t.l,t)J/Aij‘ 


(22a) 
( 22b) 


(A  Gr)  j 


•9 

C  A*,  j  «?,.*  .*—{7,  .i.i,.,)  —  A*. f,t>1 — G,.i,t.s)]/Ap* 


i.k-G..i-i  >-Gl.,.l.ll  +  Cl.,.,.t .t)/(A£Aij) 
—  (G,.j,i,—G,, i.,..,  —G,  ).,.*  +  G,,, /(.At)  AC) 
Gi ,  ,  *+Ci.l.,.»_l)/A *  AC) 


If  we  add  the  obvious  form  G„ 


(22c) 
( 22d) 
(  2  2  e } 
(22f) 


e  At  G„ *  -eAt  <«/  G i , + vgG , , + wh  G: ,) 

A|‘q~[  ~  A*'  '  »  +  -^-(Gl.t.k-Gi,  ,.t- 

{**»!-!.» +Go,i-i;»)+  ££*  <G *,,t'l>-G,.,,k-G*i,t'k.i+G,.t.i.l)  j 


(23) 
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In  all  situations,  the  velocity  components  «.  v.  w  ail  use  the 
old  values  of  the  central  difference.  If  we  take  corresponding 
difference  forms  and  substitute  them  into  equations  (17)  ana  (19)  • 

In  the  direction-,  we  obtain  the  three  diagonal  algebraic  equation 
set  . 

yi.G:. /?,<v; JM  '  <2\) 

We  use  follow  up  methods  to  solve  this. 


V.  The  Treatment  of  Boundary  Conditions 


The  object  surface  condition  (8)  for  the  upper  surface  becomes 
( CO*  «  h{°dx  cl)  '/.  d$V  C?:,-0  (25 


In  the  interior  of  each  of  the  wing  cross  sections,  we  introduce  the 
virtual  grid  points  <i ,j,k)  (See  Pig.  2). 


Pig.  2  The  Relationship  Between  the  Porm  of  the  Wing  Cross  Section 
and  the  Grid  Network  (Upper  Surface)  1 .  Upper  Surface 
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If  we  take 
virtual  points 


<M».  (r„;,  <7 


and  make  a  Taylor  expansion  on  the 


fi>ii-irffi.i,i.i,i  .  Oh  —  TJi-i)  in  n  in 


C„-  -?£«■• izi  lZj^‘,iA±^L f  (Jn _ ,  j  J*iimjLZlGu^±GJ.JrJy 

.  2Ajj  Ag, 

*  \ 

(7,.a  ***^<.  l-i.  *-l  .1  V*  *7:- 1  //;  p  ip  ,  p»  i 

f*  2A£  ^  2A^4{ 


(26) 


If  we  then  take  the  expansion  form  (26)  and  substitute  it  into  (25), 
it  is  possible  to  obtain  an  expression  for  the  perturbation  velocity  <7 
potential  at  the  virtual  points.  If  flow  movements  at  the  points  O’,  j,  k) 
are  of  supersonic  speed,  and,  *>o,  f  then,  the 

value  of  C..,-,.*  is  still  required  in  the  difference 

formula.  At  this  time,  **<•»-»•»  "Is  determined  from  linear 

extrapolation  to  be 


(27) 


The  calculation  of  the  perturbation  velocity  potential  at  the  virtual 
points  also  involves  the  introduction  of  a  relaxation  process.  First 
of  all,  we  make  use  of  the  last  iterative  substitution  value  to  139 

calculate  the  old  value  of  •  After  that,  following 

the  execution  of  a  series  of  relaxation  solutions,  we  make  use  of  as 
many  new  values  as  possible  and  recalculate  new  values  for  ^ 


Boundary  conditions  for  the  lower  surface  can  be  handled  in  tne 
same  way 

On  the  Trefftz  plane,  in  the  coordinate  system  (£,n,C) 

the  Laplace  equation  becomes 

0(06.), +*</»(?:>;-  0 


(23) 
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In  areas  of  wake  turbulence,  the  perturbation  velocity  potential  •<?  is 
not  continuous, 

G<!,+0,£>-G(*.-0,£>«r(C)  (29) 

..VI .  Result  Analysis 

We  first  make  use  of  a  rectangular  wing  and  two  sweptback  wings 
(NACA'  64A  010  and  M6X-to  carry  out  transonic  linear  flow  calculations. 
Pig.  3  and  Pig.  4  show  that  the  results  of  calculations  for  the  two 
types  of  configurations  of  the  NACA  0012  rectangular  wing  (  (M 0.75,  a=*2e 
and  JI/.-0.85,  )  when  compared  to  the  results  of 

James on  calculations,  show  precise  matching  of  the  position  of  the 
shock  waves  and  a  basic  correspondence  in  pressure  distribution. 


Pig.  3  Comparison  of  Pressure  Distributions  for  NACA  001 2  Rectangular 
Wing  (  '*  =6)  (  M~»o.75,«-r  )  i.  This  Article  2.  Upper  Surface 

3.  Lower  Surface  4.  Calculation  Results 

Pig.  5  shows  that  the  NACA  64A  010  equal  chord  back  swept  wing 
(forward  edge  back  sweep  angle  30  degrees,  expansion  chord  ratio  4.8). 
When  the  results  of  calculations  under  the  conditions  0.9,a~ie 

and  the  results  of  calculations  with  finite  body  integration 
methods  as  in  reference  [4]  and  Jameson  are  compared  to  each  other, 
pressure  distributions  and  shock  wave  locations  are  both  in  relatively 
good  agreement. 
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Pig.  4  Pressure  Distribution  for  NACAD012  Rectangular  Wing  (  <* 

(  (MM-o.js.a.o'j  )  i  .  Expansion  Chord  Position  2.  Basic 

Calculation  3*  Experiment 

Pig.  5  Pressure  Distribution  on  Surface  of  NACA  64  A  010  Back  Swept 
Wing  (  a  =  4.8)  (  »•.»,••!•>  )  y  .  Reference  4  2.  Basic 

Calculation 

150 

Finally,  we  calculated  the  ONERAM  6  wing  in  the  M  •-  3.06' 

configuration.  The  front  edge  hack  sweep  angle  of  the  wing  is  30 
degrees.  The  taper  ratio  is  0.5624.  The  expanision  chord  ratio  is 
3.8.  The  solution  for  this  has  a  double  shock  wave.  It  contains 
oblique  shock  waves  of  from  transonic  speeds  to  supersonic  speeds 
following  a  main  shock  wave  which  is  almost  perpendicular  to  the 
oncoming  flow.  The  two  shock  waves,  in  the  vicinity  of  the  wing  tips, 
merge  into  a  single  shock  wave. 

Fig.  6  shows  the  pressure  coefficient  distributions  for  four 
different  wing  cross  sections  (0.1 63 »  0.47,  0.614,  0.933).  Moreover, 
it  carries  out  a  comparison  for  the  calculated  results-both 
experimental  and  Jameson-for  adjacent  cross  sections.  Looking  at  it 
as  a  whole,  these  results  are  fairly  good.  Except  for  the  area  of  the 
forward  edge  on  the  top  surface,  there  is  entire  agreement  in  pressure 
distributions  and  experimentation  at  other  places.  After  an 
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adjustment  in  the  second  shock  wave  through  comparison  with  the 
experimental  results,  we  consider  that  there  is  no  additional 
viscosity.  Therefore,  the  second  shock  waves-hefore  adjustment  by 
comparison  of  the  basic  results  with  Jameson-are  even  more  reasonable. 

However,  due  to  the  fact  that  the  grid  which  is  employed  in  these 
calculations  is  relatively  coarse  (72x32x24)-it  is  particularly  true 
that  the  grid  points  in  the  chord  direction  are  rare-there  is  a 
flattening  effect  on  the  shock  waves.  In  particular,  the  influence  on 
the  first  shock  wave  passing  from  transonic  to  supersonic  speeds  is 
relatively  large.  This  causes  results  which  are  dissimilar  from 
Jameson  (His  grid  was  192x24x32.),  and,  in  the  area  of  the  forward 
edge,  there  is  even  better  agreement  with  experimental  values.  On  the 
basis  of  the  currently  existing  calculation  results  for  the  two 
grids- the  coarse  (36x16x12)  and  the  fine  (72x32x24) -we  estimate  that  a 
continued  increase  in  the  density  of  grid  points  in  the  chord 
direction  will  cause  the  pressure  distribution  and  experimental 
results  in  the  vicinity  of  the  forward  edges  to  agree  very  wsl" 


Fig.  6  ONER AM  6  Wing  Pressure  Coefficient  Distribution  Compa~' 
)  1.  Experiment  2.  Original  Calculation 
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The  cross  section  surface  of  the  root  portion  is  difficult  to 
calculate  accurately.  The  Jameson  method  has  made  considerable 
improvement  in  this  regard.  Moreover,  we  have  not  yet  seen  its 
results  demonstrated.  The  small  perturbation  method  which  was 
improved  in  reference  [2]  has  yielded  results.  However,  due  to  the 
fact  that  the  equations  which  are  employed  in  its  use  are  not  the 
same,  it  is  only  possible  to  see  from  Pig-  7  the  coincidence  in'  trend 
of  the  two.  The  Euler  equations  of  Jameson  solve  for  the  pressure 
distributions  at  99%  of  the  cross  section  locations  along  the  length 
of  the  development  [3]  (reference  number  unclear).  At  the  same  time, 
this  reference  provides  experimental  results.  A  comparison  of  the 
original  calculation  results  and  these  results  is  shown  in  Pig.  8.  It 
is  obvious  that  the  original  results  are  even  closer  to  experimental 
values . 

The  calculations  demonstrate  that,  in  calculations  done  in‘ 
relatively  dense  grid  systems,  the  mutual  ratio  between  increment 
lengths  A.x.  Ay  and  the  selection  of  the  artificial  viscosity 
coefficient  e  have  relatively  large  influences  on  the  convergence 
characteristics  and  stability  of  the  calculations.  These  calculations 
make  use  of  adjustable  values  for  e  and  the  relaxation  constant  », 
in  order  to  increase  the  reliability  of  the  calculations. 


Pig.  7  Pressure  Distribution  on  the  Symmetrical  Surface  of  the  ho 
Wing  (  ».**.«- }  i  .  Original  Calculations  2.  Small 

Perturbation 

Pig.  8  Pressure  Distribution  on  M6  Wing  Cross  Section  in  the  0.99 
Development  Direction  (  ,M  m  )  i  ,  Original 

Calculation  3*  Experiment  4*  Euler  Equation  Solution 

. .  •:  .  .  ,  ......  ...  .  .  ....  99  . 


VII.  Conclusion 


.The  calculations  which  we  currently  have  demonstrate  that  the 
method  for  handling  the  boundary  conditions  on  the  object  surface  of 
the  flow  around  the  Carlson  two-dimensional  airfoil  is  suitable,  in' 
principle,  for  the  calculation  of  three-dimensional  airfoils.  At  the 
same  time,  they  also  demonstrate  that,  in  orthogonal  coordinate 
systems,  it  is  possible  to  obtain  relatively  accurate 
three-dimensional  transonic  flow  solutions. 
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NUMERICAL  CALCULATIONS  OF  INVISCip  TRANSONIC 

FLOWS  OVER  WINGS 

Chen  Zuobin  Yao  Furu  Zhang  Yulun 
(China  Aerodynamic  Research  and  Development  Center) 

Abstract 

A  numerical  method  for  calculating  steady  inviscid  transonic  flows 
over  wings  is  introduced  in  this  paper.  The  fully  potential  equation  is 
chosen  as  the  mathematical  model.  The  infinite  physical  domain  is  trans¬ 
formed  to  the  finite  computing  domain  and  the  swept  wing  is  transformed 
to  the  rectangle  wing,  using  proper  mathematical  transformation.  The 
equation  is  discretized  in  computing  domain  using  mixed  finite  scheme 
and  the  line  relaxation  method  is  used  for  solving  the  resulting  nonlinear 
algebraic  equation.  Successively  refining  mesh  makes  the  calculation  very 
economic. 
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THE  CALCULATION  OP  NON-STEADY  STATE 
VISCOUS  PLOWS  BY  N-S  EQUATION  AROUND 
JOUKOWSKY  AIRFOILS  HAVING  ANGLE  OP  ATTACK 


Du  Dirong  Li  Jin 


(Northwestern  Poly technical  University) 


SUMMARY  This  article  makes  use  of'N-3  equations  to  calculate 
viscous  flows  around  Joukowsky  airfoils  with  angles  of  attack  of  0 
degrees  and  42  degrees.  The  thickness  of  the  symmetrical  airfoil  is 
25 .3# .  The  non-dimensional  grid  time  intervals  are  *  0  to  f=!0. 

The  calculation  results  include:  flow  function  distribution,  wing 
surface  pressure  distribution,  and  lift  coefficients  as  well  as  drag 
coefficients. 

I.  Introduction 

When  aircraft  are  flying  with  large  angles  of  attack,  one  always 
has  a  considerable  number  of  vortical  flows  and  separati  on  flows 
created.  The  occurrence  of  the  related  steady  state  vortices  and  the 
development  process  as  well  as  such  questions  as  the  aerodynamic 
characteristics  of  conditions  with  separation  flows  all  require  the 
doing  of  more  advanced  research  and  discussion..  Because  of  the  rapid 
development  and  broad  applications  of  electronic  computers,  and  only 
because  of  them,  has  it  become  gradually  possible  to  do  research  on 
these  problems  and  find  solutions  to  them,  [l  ,4— 6 j 

This  article  intends  to  make  use  of  R-S  equations  on  the  flow 
configurations  when  non-compressible  gases  flow  around  airfoils  as 
well  as  on  the  processes  of  formation  and  development  of  vortices  in 
order  to  carry  out  numerical  value  calculation  research.  We  make 
direct  use  of  difference  methods  to  calculate  the  initial  movements  of 
Joukowsky  airfoils  starting  from  a  static  configuration  and  moving 
suddenly  to  constant  velocity  V  .  In  numerical  value  calculation 
examples,  the  thickness  corresponding  to  symmetrical  airfoils  is 
25*8$;  the  angles  of  attack  are  0  degrees  and  42  degrees;  and,  the 
Reynolds  number  is  20.  In  these  calculations,  the  non-dimensional 
grid  time  intervals  are  from  i---o  to  >  .  The  results  of 

calculations  include  graph  forms  for  flow  movements,  wing  surface 
pressure  distributions,  and  coefficients  for  lift  and  drag. 
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The  general  procedure  for  these  calculations  is  as  follows.  We 
first  make  use  of  a  Joukowsky  transform  to  take  the  plane,  physical 
surface  of  a  flow  movement  <*, y)  and  transform  it  to  plane  U,n) 

On  plane  (£.  ,  -the  wing  surface  "becomes  a  unit  circle.  At  the 

same,  time,  on  the  <«*  D  plane,  we  carry  out  a  rotational 
transformation,  causing  the  polar  axis  and  the  direction  of  the 
incoming  flow  to  coincide  in  polar  coordinate  system  (r'  )  after 

the  transformation.  In  order  to  make  the  calculations  simple,  on  the  (r,  0,> 
plane,  we  carry  out  an  extension  transformation:  r-.-e" 

The  purpose  of  this  transformation  is  to  make  the  uniform 
distribution  difference  grid  on  the  plane  <P*0t)  able  to  correspond 
to  the  grid  distribution,  which  is  loose  on  the  outside  and  dense  on 
the  inside,  on  the  physical  plane  .  At  this  time,  the  N-3 

equations  and  boundary  conditions  on  plane  (■*.  are  gradually 
transformed  onto  the  plane  0»,  0.)  according  to  the  principles 
described.  In  the  calculations,  we  selected  for  use  a  non-conserv.r i 
form  of  difference  equation  For  time  increment  derivative,  we 
selected  for  use  the  forward  difference.  For  the  space  derivative,  we 
chose  to  use  a  central  difference-  In  handling  the  difference  in 
boundary  conditions,  the  space  derivative  makes  use  of  a  single  side 
difference.  Intercept  error  is  0 (/«)  .  In  numerical  value 

calculations,  the  initial  flow  function  value  jpjJt  is  taken  to  be  an 
inviscid  potential  flow  solution.  Concerning  calculations  o'  the 
difference  equations  obtained  from  the  transfer  of  vorticity 
displacement  equations,  we  made  use  of  the  time  increment  driven 
method.  Concerning  difference  equations  obtained  from  Poisson 
equations,  we  made  use  of  the  superrelaxation,  iterative  substitution 
method  to  carry  out  the  calculations . 


II.  Coordinate  Transformations  and  Formulae  Derivation  194 

In  problems  relating  to  two  dimensional,  viscous  non-steady 
state,  non-compressible  flows,  the  N-S  equation  can  be  written  as 
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du  dv 
d'x  dy 


ot  Ox  dy 


iL  4.  L 

ax 


V'u 


dv 

dT 


+a~-  +w 
ox 


dv 

dy 


ty_  +  JL 

dy  +  Re 


We  select  for  use  the  flow  function  ^  and  the  vorticity  ©  to  be  the 
new  variables.  Then,  the  equations  above  become 

+a-V«*>5e 

—  «d 


We  choose  to  use  a  series  of  coordinate  transforms  to  simplify  the 
process  of  numerical  value  calculation.  Moreover,  we  make  use  of  the 
Joukowsky  transform  to  take  the  airfoil  and  change  it  into  a  unit 
circle. 

*“(|+<>  (i+ ({+0‘+uT) 


In  the  numerical  example,  we  choose  /= 0.8,  e«o.2.  At  this 

time,  on  the  plane  <••'»!/>  ,  the  corresponding  thickness  of  the 

Joukowsky  airfoil  is  25. 8£,  which  is  changed,  on  the  (£,  tp  plane, 
into  a  unit  circle. 

The  equations  discussed  above  in  which  W  and  o  vrere  taken  as 
new  unknowns  are  changed  to  be: 
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In  the  equations, 


W'V.'a 

H«vJ«y=-a 


«-(-£)'4£  )'•  41 


,  .  [d'l'  dV\ 


Moreover,  on  the  basis  of 


£  —  r  cos  iJ,  )/-•  r  tin  0 


it  is  possible  to  obtain  the  equation  form  U  polar  coordinates,  which 


is: 


-  +I/‘(u,  da  +u.  I-  dI_  +  J  J  1  d‘_\ 

df  \  dr  +  r  dtf  /  /te  \  dr*  +  r  dr^  r*  dflr /„ 

yr.(  d1  I  d  1  d’  \ 


When  the  angle  of  attack  is  not  zero,  we  select  for  use  on  this  a 
rotational  transform  o,**9-a 

Then,  the  form  of  the  equation  on  the  plane  'r>0,)  is; 

C,+«>  [».-'>  +..  -L  J“  ] 

aa  1  Hi  f  0*0  ,|  da  .  1  d'ol 

H?  ld7*'frdr+  vr  cWJ 

It  Hr  «  4  ti)  \<W  .  1  OV  .  1  d’VM 
"  '  •*'  f  >  I  Sr*'  +  r  J,  +  ,i 


Making  use  of  an  extension  transform 


r*#' 
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'fl,  9t)  plane  is 


Then,  the  form  of  the  equations  on  the 

Tr"+H  <  ■*,+  ’  lur  •  &--&■  -»j 

-A  +--~(v+Sr) 

Again,  we  bring  into  play  the  variable  transform 

Then,  the  equations  change  to  become 

.  .£  - .5. .  .*•  ) 

The  formula  for  calculating  wing  surface  pressure  coefficients  is  as 
follows: 


#  "■-'•"HI* 


The  formulas  for  calculating  the  coefficients  of  lift  and  drag  are  as 
given  below: 

CDF-  -  y '  2CP(9,)t,\  (OtUP, 

CLP-j"  2CP(ff,).r; 

CDF -J”- <0.U(, 

CLF"  ~ I!’  A  oi9')u'  <*•><". 
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In  these  formulas,  sin  a+x  cos  a,  cos  a—x  tin  a 

Here,  the  direction  of  the  *•  axis  In  the  <*».  y.)  coordinate 
system  is  parallel  to  the  direction  of  the  infinitely  distant  oncoming 

flow. 

Besides  this,  the  initial  conditions  and  boundary  conditions  are 

as  stated  below: 


/<0  Q£6£2x 
tn> 6  O£0£2* 
/>  0  O£0£2« 


0<p<oo  u— 0, 

p*  0  «»0,  w«0 

p«oo  0  — 1*  sin  0lf  on  <P  —  Q 


III. 


Numerical  Value  Calculation  Methods 


and  Procedures 


We  choose  the  place  at  which  p-=p».«  to  represent  an 
Infinitely  distant  boundary. 

In  domains  on  the  calculation  plane  <p,  o,)  , 

we  lay  out  a  rectangular  coordinate  grid  with  o<t/»-rp.„,  o<'^,<2 * 
the  grid  side  length  being  h  . 

Finally,  the  difference  equations  which  we  get  frou  the  equations 
which  we  g»!.  by  taV :  ng  ®  and  0  as  variables  are: 


*  ■  *£#  1 . .  -(- • -*?-•-)  •  ( •)- 


V  2A  / 

\ 2A  j 

*  ■  (  $  ,l*  *■  -  $ 
'(""“'iT""' 


Here 
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At  this  time,  the  difference  equation  which  corresponds  to  the  Poisson 
equation  in  the  equation  set  is 


In  this  equation,  the  relaxation  factor,  and  the  superscript  i 

represents  th«  n  th' iterative  substitution. 

In  the  division  of  the  grid  net,  the. radial  direction  is  45-  The 
circumferance  is  60.  is  selected  at  the  position  where  it  o  3.t/2 

The  exterior  bound  .> -j  separation  from  the  object  surface  has 
a  chord  length  of  approximately  33*34. 

Numerical  value  calculation  procedures  are  as  shown  below.  1  97 

In  calculating  H ;*  ,  the  given  boundary  placement  conditions  and 
wing  surface  conditions  provide  the  initial  values  i'or  <f>,,  and  a’>i  . 
Subsequently,  we  calculate  the  value  of  for  wing  surfaces. 

Again,  we  solve  for  the  value  of  &,t  ■  t  the  next  instcut,  and, 
following  that,  we  solve  for  Uc  value  of  at  that  same  next 
Instant.  On  the  flow  movem.-  11  .ua  printed  for  the  given  instant, 
we  also  calculate  the  values  for  CP,  CUP,  CLP,  CDF,  CLP 

In  order  to  obtain  tne  i^^uircd  precise  huh,  it  I  s  possible  to 
obtain,  from  the  calculated  value  of  ,  a  calculation  for  the  value 
of  w*».  on  wing  surfaces.  Again,  w  do  substitution  calculations  in 
order  to  facilitate  the  calculations  for  final  flow  diagrams  and  for 
aerodynamic  coefficient  values. 

In  the  process  of  these  calculations,  one  becomes  involved  in  the 
question  of  the  selection  of  initial  values,  the  question  of 
determining  exterior  boundary  conditions  for  the  lower  flows,  as  well 
as  being  involved  in  the  questions  of  grid  size  and  coordinate 
convergence.  Each  of  these  questions  requires  separate  consideration. 

When  doing  calculations  in  flow  problems  concerning  the  patterns 
of  viscous,  non-compressible  fluids  as  they  pass  objects  after  having 
been  suddenly  launched  into  motion,  we  make  use  of  flow  solutions  for 
an  ideal  fluid  (  the  amount  of  circulation  is  zero),  and  it  is 
reasonable  and  workable  to  take  these  as  initial  values  (because,  at 
this  time,  the  influence  of  viscosity  has  still  not  had  an  effect). 

On  this  foundation,  the  calculations  as  a  whole  go  smoothly  and  well. 
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Because  of  this,  to  the  extent  that  we  are  concerned  with  flow  pattern 
problems  in  viscous  gases  with  differing  angles  of  attack,  Lc  is 
appropriate  to  select  for  use  the  corresponding  potential  flow 
solutions  as  initial  values.  In  the  case  of  the  Joukowsky  airfoil,  . 
potential  solutions  for  angles  of  attack  make  it  easy  to  solve  for  and 
analyze  the  solutions.  As  far  as  airfoils  in  general  are  concerned, 
there  is  normally  no  conventional  or  appropriate  formula  expressing 
the  analyzed  solution.  In  this  case,  it  is  possible,  on  the 
foundation  of  the  relationship  <£ii=©,i=0  ,  to  adopt  a 

relaxation  iterative  substitution  method  to  solve  correspond ing  diff¬ 
erence  equations  for  Laplace  equations.  Moreover,  the  given 
circulation  is  zero.  In  this  way,  the  initial  values  which  are  solved 
for  are  actually  nothing  more  than  numerical  value  approximations  of 
potential  flow  solutions.  Of  course,  when  the  angle  of  attack  is 
zero,  the  problem  is  relatively  simpler.  This  is  particularly  true  of 
symmetrical  airfoils.  Prom  the  symmetrical  characteristics  of  the 
flow  movements  on  the  upper  and  lower  surfaces,  it  is  possible  to  make 
a  solution.  However,  in  a  situation  where  the  angle  of  attack  for  the 
airfoil  is  set  at  zero,  is  usually  also  necessary  to  make  use  of  the 
condition  that  the  circulation  is  zero  in  order  to  calculate  initial 
values. 

Here,  we  take  the  wing  surface  boundary  and  designate  it  as  the 
interior  boundary  of  the  calculation  domain. 

This  closed  curve  is  designated  as  the  exterior  p-pw.  or^<2.T 
boundary  of  the  domain  of  the  calculations.  Concerning  the 
determination  of  the  exterior  boundary  conditions  in  the  lower  flows, 
different  theorists  each  have  their  own  methods  for  doing  the 
determination  (see  references  [2]  and  [3])*  We  recognize  that,  in 
the  situation  in  which  the  exterior  boundaries  of  the  domain  of 
calculation  are  chosen  relatively  far  away  from  wing  surfaces,  it  is 
possible  to  make  initial  use  of  the  free  flow  conditions  as  the 
exterior  boundaries  of  the  lower  flows  in  order  to  carry  out  trie 
calculations.  Because  of  this,  at  this  time,  wake  vortices  have  still 
not  dragged  out  the  exterior  boundaries.  When  t  reaches  a  specified 
numerical  value  of  a  certain  size,  then,  wake  vortices  have  already 


peached  the  vicinity  of  the  "boundaries.  .  At  just,  the  time  that  that 
happens,  we  consider  opting  for  the  use  of  relatively  precise 
differential  equation  forms  in  order  to  specify  the  conditions  for  the 
exterior  "boundaries-  of  the  lower  flows.  In  these  calculations,  one - 
makes  use  of  the  conditions  for  the  free  incoming  flow  to  be  the 
exterior  boundary  conditions  for  the  lower  flows.  The  duration  of  the 
calculations  only  reaches  to  <-io.  The  perturbation  on  the  airfoil 
from  the  influence  of  the  incoming  flow- has  still  not  reached-in  the 
form  of  vortices-the  exterior  boundaries  of  the  lower  flows. 
Simplifying  the  calculations  in  this  way,  one  sees  the  calculations 
for  problems  becoming  easy  and  their  solutions  coming  in  a  timely 
manner . 

The  selection  of  the  size  of  the  grid  network  is  not  only 
effected  by  the  influence  of  the  values  of  .  It  is  also  related  to 
the  thickness  of  the  airfoil.  In  order  to  make  an  appropriate 
selection  of  the  level  of  density  of  the- gradient ,  one  should  base  it 
on  the  extent  of  the  space  invloved  in  the  changes  which  occur  in  the 
flow  movements. 

This  article  makes  use  of  the  standard  of  absolute  convergence, 
that  is,  when  the  absolute  value  of  the  difference  between  two 
adjacent  iterative  substitution  values  of  on  grid  points  is 

I £0.0001  ,then,  the  iterative  substitution  stops. 

Because  of  the  fact  that,  when  one  is  doing  a  numerical  value 
treatment  of  boundary  conditions,  the  Intercept  error  Is  0(A),,  •  i ", 
is  the  grid  length,  its  value  in  all  of  the  calculations  is  /i*.-r/30>  o.i. 

In  cases  where  this  type  of  convergence  standard  is 
employed,  in  calculations  in  which  . /?«■*  20  f  the  number-  of  iterative 
substitutions  on  the  average  in  calculations  of  +,i  vt  each  instant 
is  approximately  20.  The  effects  of  error  on  the  r- >ults  of 
calculation  are  small,  and  calculations  can  be  carried  out  with  ease. 
When  t  is  increased,  there  is  no  clear  increase  in  the  number  of 
tte-».';ive  substitutions  for  each  instant. 


Ill 


IV.  Calculation  Results  and  Discussion 


This  article  has  made  calculations  involving  a  symmetrical 
Joukowsky  airfoil  with  thickness  25.3$  and  has  dealt  with  problems  of 
flow-forms  around  objects -in  viscous  gases  when  AV-20  <=*nd  with 
angles  of  attack  which  were  respectively  0  degrees  and  42  degrees. 

The  airfoil,  at  time  *”0 ,  is  suddenly  launched  into  motion  from  a 
dead  stop.  Afcer  that-,  it  goes  into  straight  line  movement  at  uniform 
speed.  Relatively  speaking,  we  have  done  nothing  more  than  calculate 
problems  in  the  flow  forms  of  non-steady  state,  viscous  gases  as  they 
flow  around  airfoils.  The  results  of  the  calculations  have  been 
printed  out  by  computer.  These  flow  diagrams  are  given  in  Pig,  1 
through  Pig.  4.  The  pressure  distribution  curves  are  as  given  in  Pig. 
5.  The  lift  coefficient  and  drag  coefficient  curves  are  drawn  out  on 
the  basis  of  calculation  results  and  are  given  in  Pig.  6  and  Pig.  7 
respectively.  Due  to  the  fact  that,  at  present,  it  is  <1  i  P f  i r.  to 


Pig.  1  Plow  Diagram  t’.ij 

Pig.  2  Plow  Diagram  r.« 

find  direct  comparisons  to  apply  to  experimental  results  or  numerical 
value  calculation  results  which  come  from  differing  situations.  It 
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Fig.  3  Flow  Diagram  «-4i*.  r-7 

Fig.  4  Flow  Diagr \  >  Re-2o.  «-«•,  t-io 

is  only  with  the  viscous  flow  diagrams  for  flows  around  long  and 
slender  elliptical  columns  in  reference  [2]  that  we  carried  out 
comparisons.  These  demonstrated  that  the  calculations  were  basical-!, 
reasonable. 


Fig.  5  Pressure  Distribution  *«•**,«•  < *\r«r  •  cr  i  o,./  •'  <i' 

Surface)  OCP  (Lo*sr  Wing  Surface;  *«•*•. 

Fig.  6  Lift  Coefficients 


Pig. 


7  Drag  Coefficients  K* *»,«*«• 
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A  NUMERICAL  METHOD  OF  UNSTEADY  FLOW 
FOR  SOLVING  NAVIER-STOKES  EQUATION 

- PAST  A  JOUKOWSKY  AIRFOIL 

WITH  ANGLE  OF  ATTACK 


Du  Dirong  Li  Jin 
(fforthwcitcra  Polytechnics!  University) 

Abstract 

The  present  paper  is  the  numerical  solution  of  the  unsteady  incompres¬ 
sible  viscous  flow  over  Joukowsky  airfoil  by  Navier-stokes  equation  with 
Mgl*  *f  attack  of  0*  and  42*.  The  thickness  ratio  of  the  symmetrical 
uirfoil  is  2S,R*».  Nondimcnsional  time  interval  is  from  t  —  G  to  t**10.  Nu¬ 
merical  results  of  interest,  such  as  stream  function  distribution,  surface 
pressure  distribution,  lift  uud  drag  coefficients  arc  computed. 
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